1 2 #include <private/pcimpl.h> /*I "petscpc.h" I*/ 3 #include <petscdmcomposite.h> /*I "petscdmcomposite.h" I*/ 4 #include <petscdmcomplex.h> 5 6 const char *const PCFieldSplitSchurPreTypes[] = {"SELF","DIAG","USER","PCFieldSplitSchurPreType","PC_FIELDSPLIT_SCHUR_PRE_",0}; 7 const char *const PCFieldSplitSchurFactorizationTypes[] = {"DIAG","LOWER","UPPER","FULL","PCFieldSplitSchurFactorizationType","PC_FIELDSPLIT_SCHUR_FACTORIZATION_",0}; 8 9 typedef enum { 10 PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG, 11 PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER, 12 PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER, 13 PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL 14 } PCFieldSplitSchurFactorizationType; 15 16 typedef struct _PC_FieldSplitLink *PC_FieldSplitLink; 17 struct _PC_FieldSplitLink { 18 KSP ksp; 19 Vec x,y; 20 char *splitname; 21 PetscInt nfields; 22 PetscInt *fields; 23 VecScatter sctx; 24 IS is; 25 PC_FieldSplitLink next,previous; 26 }; 27 28 typedef struct { 29 PCCompositeType type; 30 PetscBool defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */ 31 PetscBool splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */ 32 PetscBool realdiagonal; /* Flag to use the diagonal blocks of mat preconditioned by pmat, instead of just pmat */ 33 PetscInt bs; /* Block size for IS and Mat structures */ 34 PetscInt nsplits; /* Number of field divisions defined */ 35 Vec *x,*y,w1,w2; 36 Mat *mat; /* The diagonal block for each split */ 37 Mat *pmat; /* The preconditioning diagonal block for each split */ 38 Mat *Afield; /* The rows of the matrix associated with each split */ 39 PetscBool issetup; 40 /* Only used when Schur complement preconditioning is used */ 41 Mat B; /* The (0,1) block */ 42 Mat C; /* The (1,0) block */ 43 Mat schur; /* The Schur complement S = A11 - A10 A00^{-1} A01 */ 44 Mat schur_user; /* User-provided preconditioning matrix for the Schur complement */ 45 PCFieldSplitSchurPreType schurpre; /* Determines which preconditioning matrix is used for the Schur complement */ 46 PCFieldSplitSchurFactorizationType schurfactorization; 47 KSP kspschur; /* The solver for S */ 48 PC_FieldSplitLink head; 49 PetscBool reset; /* indicates PCReset() has been last called on this object, hack */ 50 PetscBool suboptionsset; /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */ 51 } PC_FieldSplit; 52 53 /* 54 Notes: there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of 55 inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the 56 PC you could change this. 57 */ 58 59 /* This helper is so that setting a user-provided preconditioning matrix is orthogonal to choosing to use it. This way the 60 * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */ 61 static Mat FieldSplitSchurPre(PC_FieldSplit *jac) 62 { 63 switch (jac->schurpre) { 64 case PC_FIELDSPLIT_SCHUR_PRE_SELF: return jac->schur; 65 case PC_FIELDSPLIT_SCHUR_PRE_DIAG: return jac->pmat[1]; 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 #undef __FUNCT__ 74 #define __FUNCT__ "PCView_FieldSplit" 75 static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer) 76 { 77 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 78 PetscErrorCode ierr; 79 PetscBool iascii; 80 PetscInt i,j; 81 PC_FieldSplitLink ilink = jac->head; 82 83 PetscFunctionBegin; 84 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 85 if (iascii) { 86 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with %s composition: total splits = %D, blocksize = %D\n",PCCompositeTypes[jac->type],jac->nsplits,jac->bs);CHKERRQ(ierr); 87 ierr = PetscViewerASCIIPrintf(viewer," Solver info for each split is in the following KSP objects:\n");CHKERRQ(ierr); 88 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 89 for (i=0; i<jac->nsplits; i++) { 90 if (ilink->fields) { 91 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 92 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 93 for (j=0; j<ilink->nfields; j++) { 94 if (j > 0) { 95 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 96 } 97 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 98 } 99 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 100 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 101 } else { 102 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 103 } 104 ierr = KSPView(ilink->ksp,viewer);CHKERRQ(ierr); 105 ilink = ilink->next; 106 } 107 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 108 } else { 109 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 110 } 111 PetscFunctionReturn(0); 112 } 113 114 #undef __FUNCT__ 115 #define __FUNCT__ "PCView_FieldSplit_Schur" 116 static PetscErrorCode PCView_FieldSplit_Schur(PC pc,PetscViewer viewer) 117 { 118 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 119 PetscErrorCode ierr; 120 PetscBool iascii; 121 PetscInt i,j; 122 PC_FieldSplitLink ilink = jac->head; 123 KSP ksp; 124 125 PetscFunctionBegin; 126 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 127 if (iascii) { 128 ierr = PetscViewerASCIIPrintf(viewer," FieldSplit with Schur preconditioner, blocksize = %D, factorization %s\n",jac->bs,PCFieldSplitSchurFactorizationTypes[jac->schurfactorization]);CHKERRQ(ierr); 129 ierr = PetscViewerASCIIPrintf(viewer," Split info:\n");CHKERRQ(ierr); 130 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 131 for (i=0; i<jac->nsplits; i++) { 132 if (ilink->fields) { 133 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 135 for (j=0; j<ilink->nfields; j++) { 136 if (j > 0) { 137 ierr = PetscViewerASCIIPrintf(viewer,",");CHKERRQ(ierr); 138 } 139 ierr = PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);CHKERRQ(ierr); 140 } 141 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 142 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 143 } else { 144 ierr = PetscViewerASCIIPrintf(viewer,"Split number %D Defined by IS\n",i);CHKERRQ(ierr); 145 } 146 ilink = ilink->next; 147 } 148 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for A00 block \n");CHKERRQ(ierr); 149 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 150 if (jac->schur) { 151 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 152 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 153 } else { 154 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 155 } 156 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 157 ierr = PetscViewerASCIIPrintf(viewer,"KSP solver for S = A11 - A10 inv(A00) A01 \n");CHKERRQ(ierr); 158 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 159 if (jac->kspschur) { 160 ierr = KSPView(jac->kspschur,viewer);CHKERRQ(ierr); 161 } else { 162 ierr = PetscViewerASCIIPrintf(viewer," not yet available\n");CHKERRQ(ierr); 163 } 164 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 165 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 166 } else { 167 SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name); 168 } 169 PetscFunctionReturn(0); 170 } 171 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCFieldSplitSetRuntimeSplits_Private" 174 /* Precondition: jac->bs is set to a meaningful value */ 175 static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc) 176 { 177 PetscErrorCode ierr; 178 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 179 PetscInt i,nfields,*ifields; 180 PetscBool flg; 181 char optionname[128],splitname[8]; 182 183 PetscFunctionBegin; 184 ierr = PetscMalloc(jac->bs*sizeof(PetscInt),&ifields);CHKERRQ(ierr); 185 for (i=0,flg=PETSC_TRUE; ; i++) { 186 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 187 ierr = PetscSNPrintf(optionname,sizeof optionname,"-pc_fieldsplit_%D_fields",i);CHKERRQ(ierr); 188 nfields = jac->bs; 189 ierr = PetscOptionsGetIntArray(((PetscObject)pc)->prefix,optionname,ifields,&nfields,&flg);CHKERRQ(ierr); 190 if (!flg) break; 191 if (!nfields) SETERRQ(((PetscObject) pc)->comm,PETSC_ERR_USER,"Cannot list zero fields"); 192 ierr = PCFieldSplitSetFields(pc,splitname,nfields,ifields);CHKERRQ(ierr); 193 } 194 if (i > 0) { 195 /* Makes command-line setting of splits take precedence over setting them in code. 196 Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would 197 create new splits, which would probably not be what the user wanted. */ 198 jac->splitdefined = PETSC_TRUE; 199 } 200 ierr = PetscFree(ifields);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "PCFieldSplitSetDefaults" 206 static PetscErrorCode PCFieldSplitSetDefaults(PC pc) 207 { 208 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 209 PetscErrorCode ierr; 210 PC_FieldSplitLink ilink = jac->head; 211 PetscBool flg = PETSC_FALSE,stokes = PETSC_FALSE; 212 PetscInt i; 213 214 PetscFunctionBegin; 215 if (!ilink) { 216 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 217 if (pc->dm && !stokes) { 218 PetscBool dmcomposite, dmcomplex; 219 ierr = PetscTypeCompare((PetscObject)pc->dm,DMCOMPOSITE,&dmcomposite);CHKERRQ(ierr); 220 ierr = PetscTypeCompare((PetscObject)pc->dm,DMCOMPLEX,&dmcomplex);CHKERRQ(ierr); 221 if (dmcomposite) { 222 PetscInt nDM; 223 IS *fields; 224 DM *dms; 225 ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 226 ierr = DMCompositeGetNumberDM(pc->dm,&nDM);CHKERRQ(ierr); 227 ierr = DMCompositeGetGlobalISs(pc->dm,&fields);CHKERRQ(ierr); 228 ierr = PetscMalloc(nDM*sizeof(DM),&dms);CHKERRQ(ierr); 229 ierr = DMCompositeGetEntriesArray(pc->dm,dms);CHKERRQ(ierr); 230 for (i=0; i<nDM; i++) { 231 char buf[256]; 232 const char *splitname; 233 /* Split naming precedence: object name, prefix, number */ 234 splitname = ((PetscObject)pc->dm)->name; 235 if (!splitname) { 236 ierr = PetscObjectGetOptionsPrefix((PetscObject)dms[i],&splitname);CHKERRQ(ierr); 237 if (splitname) { 238 size_t len; 239 ierr = PetscStrncpy(buf,splitname,sizeof buf);CHKERRQ(ierr); 240 buf[sizeof buf - 1] = 0; 241 ierr = PetscStrlen(buf,&len);CHKERRQ(ierr); 242 if (buf[len-1] == '_') buf[len-1] = 0; /* Remove trailing underscore if it was used */ 243 splitname = buf; 244 } 245 } 246 if (!splitname) { 247 ierr = PetscSNPrintf(buf,sizeof buf,"%D",i);CHKERRQ(ierr); 248 splitname = buf; 249 } 250 ierr = PCFieldSplitSetIS(pc,splitname,fields[i]);CHKERRQ(ierr); 251 ierr = ISDestroy(&fields[i]);CHKERRQ(ierr); 252 } 253 ierr = PetscFree(fields);CHKERRQ(ierr); 254 for (ilink=jac->head,i=0; ilink; ilink=ilink->next,i++) { 255 ierr = KSPSetDM(ilink->ksp,dms[i]);CHKERRQ(ierr); 256 ierr = KSPSetDMActive(ilink->ksp,PETSC_FALSE);CHKERRQ(ierr); 257 } 258 ierr = PetscFree(dms);CHKERRQ(ierr); 259 } else if (dmcomplex) { 260 /* Define fields */ 261 PetscSection section, sectionGlobal; 262 PetscInt *fieldSizes, **fieldIndices; 263 PetscInt nF, f, pStart, pEnd, p; 264 265 ierr = DMComplexGetDefaultSection(pc->dm, §ion);CHKERRQ(ierr); 266 ierr = DMComplexGetDefaultGlobalSection(pc->dm, §ionGlobal);CHKERRQ(ierr); 267 ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); 268 ierr = PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt *,&fieldIndices);CHKERRQ(ierr); 269 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 270 for(f = 0; f < nF; ++f) { 271 fieldSizes[f] = 0; 272 } 273 for(p = pStart; p < pEnd; ++p) { 274 PetscInt gdof; 275 276 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 277 if (gdof > 0) { 278 for(f = 0; f < nF; ++f) { 279 PetscInt fcdof; 280 281 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); 282 fieldSizes[f] += fcdof; 283 } 284 } 285 } 286 for(f = 0; f < nF; ++f) { 287 ierr = PetscMalloc(fieldSizes[f] * sizeof(PetscInt), &fieldIndices[f]);CHKERRQ(ierr); 288 fieldSizes[f] = 0; 289 } 290 for(p = pStart; p < pEnd; ++p) { 291 PetscInt gdof, goff; 292 293 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 294 if (gdof > 0) { 295 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); 296 for(f = 0; f < nF; ++f) { 297 PetscInt fcdof, fc; 298 299 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); 300 for(fc = 0; fc < fcdof; ++fc, ++fieldSizes[f]) { 301 fieldIndices[f][fieldSizes[f]] = goff++; 302 } 303 } 304 } 305 } 306 for(f = 0; f < nF; ++f) { 307 IS field; 308 const char *fieldname; 309 310 ierr = PetscSectionGetFieldName(section, f, &fieldname);CHKERRQ(ierr); 311 ierr = ISCreateGeneral(((PetscObject) pc)->comm, fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &field);CHKERRQ(ierr); 312 ierr = PetscPrintf(((PetscObject) pc)->comm, "Field %s\n", fieldname);CHKERRQ(ierr); 313 ierr = ISView(field, PETSC_NULL);CHKERRQ(ierr); 314 ierr = PCFieldSplitSetIS(pc, fieldname, field);CHKERRQ(ierr); 315 ierr = ISDestroy(&field);CHKERRQ(ierr); 316 } 317 ierr = PetscFree2(fieldSizes,fieldIndices);CHKERRQ(ierr); 318 } 319 } else { 320 if (jac->bs <= 0) { 321 if (pc->pmat) { 322 ierr = MatGetBlockSize(pc->pmat,&jac->bs);CHKERRQ(ierr); 323 } else { 324 jac->bs = 1; 325 } 326 } 327 328 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 329 if (stokes) { 330 IS zerodiags,rest; 331 PetscInt nmin,nmax; 332 333 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 334 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 335 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 336 ierr = PCFieldSplitSetIS(pc,"0",rest);CHKERRQ(ierr); 337 ierr = PCFieldSplitSetIS(pc,"1",zerodiags);CHKERRQ(ierr); 338 ierr = ISDestroy(&zerodiags);CHKERRQ(ierr); 339 ierr = ISDestroy(&rest);CHKERRQ(ierr); 340 } else { 341 if (!flg) { 342 /* Allow user to set fields from command line, if bs was known at the time of PCSetFromOptions_FieldSplit() 343 then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */ 344 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 345 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 346 } 347 if (flg || !jac->splitdefined) { 348 ierr = PetscInfo(pc,"Using default splitting of fields\n");CHKERRQ(ierr); 349 for (i=0; i<jac->bs; i++) { 350 char splitname[8]; 351 ierr = PetscSNPrintf(splitname,sizeof splitname,"%D",i);CHKERRQ(ierr); 352 ierr = PCFieldSplitSetFields(pc,splitname,1,&i);CHKERRQ(ierr); 353 } 354 jac->defaultsplit = PETSC_TRUE; 355 } 356 } 357 } 358 } else if (jac->nsplits == 1) { 359 if (ilink->is) { 360 IS is2; 361 PetscInt nmin,nmax; 362 363 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 364 ierr = ISComplement(ilink->is,nmin,nmax,&is2);CHKERRQ(ierr); 365 ierr = PCFieldSplitSetIS(pc,"1",is2);CHKERRQ(ierr); 366 ierr = ISDestroy(&is2);CHKERRQ(ierr); 367 } else SETERRQ(((PetscObject) pc)->comm,PETSC_ERR_SUP,"Must provide at least two sets of fields to PCFieldSplit()"); 368 } else if (jac->reset) { 369 /* PCReset() has been called on this PC, ilink exists but all IS and Vec data structures in it must be rebuilt 370 This is basically the !ilink portion of code above copied from above and the allocation of the ilinks removed 371 since they already exist. This should be totally rewritten */ 372 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 373 if (pc->dm && !stokes) { 374 PetscBool dmcomposite; 375 ierr = PetscTypeCompare((PetscObject)pc->dm,DMCOMPOSITE,&dmcomposite);CHKERRQ(ierr); 376 if (dmcomposite) { 377 PetscInt nDM; 378 IS *fields; 379 DM *dms; 380 ierr = PetscInfo(pc,"Setting up physics based fieldsplit preconditioner using the embedded DM\n");CHKERRQ(ierr); 381 ierr = DMCompositeGetNumberDM(pc->dm,&nDM);CHKERRQ(ierr); 382 ierr = DMCompositeGetGlobalISs(pc->dm,&fields);CHKERRQ(ierr); 383 ierr = PetscMalloc(nDM*sizeof(DM),&dms);CHKERRQ(ierr); 384 ierr = DMCompositeGetEntriesArray(pc->dm,dms);CHKERRQ(ierr); 385 for (i=0; i<nDM; i++) { 386 ierr = KSPSetDM(ilink->ksp,dms[i]);CHKERRQ(ierr); 387 ierr = KSPSetDMActive(ilink->ksp,PETSC_FALSE);CHKERRQ(ierr); 388 ilink->is = fields[i]; 389 ilink = ilink->next; 390 } 391 ierr = PetscFree(fields);CHKERRQ(ierr); 392 ierr = PetscFree(dms);CHKERRQ(ierr); 393 } 394 } else { 395 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);CHKERRQ(ierr); 396 if (stokes) { 397 IS zerodiags,rest; 398 PetscInt nmin,nmax; 399 400 ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr); 401 ierr = MatFindZeroDiagonals(pc->mat,&zerodiags);CHKERRQ(ierr); 402 ierr = ISComplement(zerodiags,nmin,nmax,&rest);CHKERRQ(ierr); 403 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 404 ierr = ISDestroy(&ilink->next->is);CHKERRQ(ierr); 405 ilink->is = rest; 406 ilink->next->is = zerodiags; 407 } else SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cases not yet handled when PCReset() was used"); 408 } 409 } 410 411 if (jac->nsplits < 2) SETERRQ1(((PetscObject) pc)->comm,PETSC_ERR_PLIB,"Unhandled case, must have at least two fields, not %d", jac->nsplits); 412 PetscFunctionReturn(0); 413 } 414 415 #undef __FUNCT__ 416 #define __FUNCT__ "PCSetUp_FieldSplit" 417 static PetscErrorCode PCSetUp_FieldSplit(PC pc) 418 { 419 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 420 PetscErrorCode ierr; 421 PC_FieldSplitLink ilink; 422 PetscInt i,nsplit,ccsize; 423 MatStructure flag = pc->flag; 424 PetscBool sorted; 425 426 PetscFunctionBegin; 427 ierr = PCFieldSplitSetDefaults(pc);CHKERRQ(ierr); 428 nsplit = jac->nsplits; 429 ilink = jac->head; 430 431 /* get the matrices for each split */ 432 if (!jac->issetup) { 433 PetscInt rstart,rend,nslots,bs; 434 435 jac->issetup = PETSC_TRUE; 436 437 /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */ 438 bs = jac->bs; 439 ierr = MatGetOwnershipRange(pc->pmat,&rstart,&rend);CHKERRQ(ierr); 440 ierr = MatGetLocalSize(pc->pmat,PETSC_NULL,&ccsize);CHKERRQ(ierr); 441 nslots = (rend - rstart)/bs; 442 for (i=0; i<nsplit; i++) { 443 if (jac->defaultsplit) { 444 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+i,nsplit,&ilink->is);CHKERRQ(ierr); 445 } else if (!ilink->is) { 446 if (ilink->nfields > 1) { 447 PetscInt *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields; 448 ierr = PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);CHKERRQ(ierr); 449 for (j=0; j<nslots; j++) { 450 for (k=0; k<nfields; k++) { 451 ii[nfields*j + k] = rstart + bs*j + fields[k]; 452 } 453 } 454 ierr = ISCreateGeneral(((PetscObject)pc)->comm,nslots*nfields,ii,PETSC_OWN_POINTER,&ilink->is);CHKERRQ(ierr); 455 } else { 456 ierr = ISCreateStride(((PetscObject)pc)->comm,nslots,rstart+ilink->fields[0],bs,&ilink->is);CHKERRQ(ierr); 457 } 458 } 459 ierr = ISSorted(ilink->is,&sorted);CHKERRQ(ierr); 460 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Fields must be sorted when creating split"); 461 ilink = ilink->next; 462 } 463 } 464 465 ilink = jac->head; 466 if (!jac->pmat) { 467 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->pmat);CHKERRQ(ierr); 468 for (i=0; i<nsplit; i++) { 469 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 470 ilink = ilink->next; 471 } 472 } else { 473 for (i=0; i<nsplit; i++) { 474 ierr = MatGetSubMatrix(pc->pmat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->pmat[i]);CHKERRQ(ierr); 475 ilink = ilink->next; 476 } 477 } 478 if (jac->realdiagonal) { 479 ilink = jac->head; 480 if (!jac->mat) { 481 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->mat);CHKERRQ(ierr); 482 for (i=0; i<nsplit; i++) { 483 ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_INITIAL_MATRIX,&jac->mat[i]);CHKERRQ(ierr); 484 ilink = ilink->next; 485 } 486 } else { 487 for (i=0; i<nsplit; i++) { 488 if (jac->mat[i]) {ierr = MatGetSubMatrix(pc->mat,ilink->is,ilink->is,MAT_REUSE_MATRIX,&jac->mat[i]);CHKERRQ(ierr);} 489 ilink = ilink->next; 490 } 491 } 492 } else { 493 jac->mat = jac->pmat; 494 } 495 496 if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR) { 497 /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */ 498 ilink = jac->head; 499 if (!jac->Afield) { 500 ierr = PetscMalloc(nsplit*sizeof(Mat),&jac->Afield);CHKERRQ(ierr); 501 for (i=0; i<nsplit; i++) { 502 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_INITIAL_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 503 ilink = ilink->next; 504 } 505 } else { 506 for (i=0; i<nsplit; i++) { 507 ierr = MatGetSubMatrix(pc->mat,ilink->is,PETSC_NULL,MAT_REUSE_MATRIX,&jac->Afield[i]);CHKERRQ(ierr); 508 ilink = ilink->next; 509 } 510 } 511 } 512 513 if (jac->type == PC_COMPOSITE_SCHUR) { 514 IS ccis; 515 PetscInt rstart,rend; 516 if (nsplit != 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"To use Schur complement preconditioner you must have exactly 2 fields"); 517 518 /* When extracting off-diagonal submatrices, we take complements from this range */ 519 ierr = MatGetOwnershipRangeColumn(pc->mat,&rstart,&rend);CHKERRQ(ierr); 520 521 /* need to handle case when one is resetting up the preconditioner */ 522 if (jac->schur) { 523 ilink = jac->head; 524 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 525 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->B);CHKERRQ(ierr); 526 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 527 ilink = ilink->next; 528 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 529 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_REUSE_MATRIX,&jac->C);CHKERRQ(ierr); 530 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 531 ierr = MatSchurComplementUpdate(jac->schur,jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->pmat[1],pc->flag);CHKERRQ(ierr); 532 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),pc->flag);CHKERRQ(ierr); 533 534 } else { 535 KSP ksp; 536 char schurprefix[256]; 537 538 /* extract the A01 and A10 matrices */ 539 ilink = jac->head; 540 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 541 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->B);CHKERRQ(ierr); 542 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 543 ilink = ilink->next; 544 ierr = ISComplement(ilink->is,rstart,rend,&ccis);CHKERRQ(ierr); 545 ierr = MatGetSubMatrix(pc->mat,ilink->is,ccis,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 546 ierr = ISDestroy(&ccis);CHKERRQ(ierr); 547 /* Use mat[0] (diagonal block of the real matrix) preconditioned by pmat[0] */ 548 ierr = MatCreateSchurComplement(jac->mat[0],jac->pmat[0],jac->B,jac->C,jac->mat[1],&jac->schur);CHKERRQ(ierr); 549 /* set tabbing and options prefix of KSP inside the MatSchur */ 550 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 551 ierr = PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,2);CHKERRQ(ierr); 552 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",jac->head->splitname);CHKERRQ(ierr); 553 ierr = KSPSetOptionsPrefix(ksp,schurprefix);CHKERRQ(ierr); 554 /* Need to call this everytime because new matrix is being created */ 555 ierr = MatSetFromOptions(jac->schur);CHKERRQ(ierr); 556 ierr = MatSetUp(jac->schur);CHKERRQ(ierr); 557 558 ierr = KSPCreate(((PetscObject)pc)->comm,&jac->kspschur);CHKERRQ(ierr); 559 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->kspschur);CHKERRQ(ierr); 560 ierr = PetscObjectIncrementTabLevel((PetscObject)jac->kspschur,(PetscObject)pc,1);CHKERRQ(ierr); 561 ierr = KSPSetOperators(jac->kspschur,jac->schur,FieldSplitSchurPre(jac),DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 562 if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) { 563 PC pc; 564 ierr = KSPGetPC(jac->kspschur,&pc);CHKERRQ(ierr); 565 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 566 /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */ 567 } 568 ierr = PetscSNPrintf(schurprefix,sizeof schurprefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 569 ierr = KSPSetOptionsPrefix(jac->kspschur,schurprefix);CHKERRQ(ierr); 570 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 571 /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */ 572 ierr = KSPSetFromOptions(jac->kspschur);CHKERRQ(ierr); 573 574 ierr = PetscMalloc2(2,Vec,&jac->x,2,Vec,&jac->y);CHKERRQ(ierr); 575 ierr = MatGetVecs(jac->pmat[0],&jac->x[0],&jac->y[0]);CHKERRQ(ierr); 576 ierr = MatGetVecs(jac->pmat[1],&jac->x[1],&jac->y[1]);CHKERRQ(ierr); 577 ilink = jac->head; 578 ilink->x = jac->x[0]; ilink->y = jac->y[0]; 579 ilink = ilink->next; 580 ilink->x = jac->x[1]; ilink->y = jac->y[1]; 581 } 582 } else { 583 /* set up the individual PCs */ 584 i = 0; 585 ilink = jac->head; 586 while (ilink) { 587 ierr = KSPSetOperators(ilink->ksp,jac->mat[i],jac->pmat[i],flag);CHKERRQ(ierr); 588 /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */ 589 if (!jac->suboptionsset) {ierr = KSPSetFromOptions(ilink->ksp);CHKERRQ(ierr);} 590 ierr = KSPSetUp(ilink->ksp);CHKERRQ(ierr); 591 i++; 592 ilink = ilink->next; 593 } 594 595 /* create work vectors for each split */ 596 if (!jac->x) { 597 ierr = PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);CHKERRQ(ierr); 598 ilink = jac->head; 599 for (i=0; i<nsplit; i++) { 600 Vec *vl,*vr; 601 602 ierr = KSPGetVecs(ilink->ksp,1,&vr,1,&vl);CHKERRQ(ierr); 603 ilink->x = *vr; 604 ilink->y = *vl; 605 ierr = PetscFree(vr);CHKERRQ(ierr); 606 ierr = PetscFree(vl);CHKERRQ(ierr); 607 jac->x[i] = ilink->x; 608 jac->y[i] = ilink->y; 609 ilink = ilink->next; 610 } 611 } 612 } 613 614 615 if (!jac->head->sctx) { 616 Vec xtmp; 617 618 /* compute scatter contexts needed by multiplicative versions and non-default splits */ 619 620 ilink = jac->head; 621 ierr = MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);CHKERRQ(ierr); 622 for (i=0; i<nsplit; i++) { 623 ierr = VecScatterCreate(xtmp,ilink->is,jac->x[i],PETSC_NULL,&ilink->sctx);CHKERRQ(ierr); 624 ilink = ilink->next; 625 } 626 ierr = VecDestroy(&xtmp);CHKERRQ(ierr); 627 } 628 jac->suboptionsset = PETSC_TRUE; 629 PetscFunctionReturn(0); 630 } 631 632 #define FieldSplitSplitSolveAdd(ilink,xx,yy) \ 633 (VecScatterBegin(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 634 VecScatterEnd(ilink->sctx,xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD) || \ 635 KSPSolve(ilink->ksp,ilink->x,ilink->y) || \ 636 VecScatterBegin(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE) || \ 637 VecScatterEnd(ilink->sctx,ilink->y,yy,ADD_VALUES,SCATTER_REVERSE)) 638 639 #undef __FUNCT__ 640 #define __FUNCT__ "PCApply_FieldSplit_Schur" 641 static PetscErrorCode PCApply_FieldSplit_Schur(PC pc,Vec x,Vec y) 642 { 643 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 644 PetscErrorCode ierr; 645 KSP ksp; 646 PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next; 647 648 PetscFunctionBegin; 649 ierr = MatSchurComplementGetKSP(jac->schur,&ksp);CHKERRQ(ierr); 650 651 switch (jac->schurfactorization) { 652 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_DIAG: 653 /* [A00 0; 0 -S], positive definite, suitable for MINRES */ 654 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 655 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 656 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 657 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 658 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 659 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 660 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 661 ierr = VecScale(ilinkD->y,-1.);CHKERRQ(ierr); 662 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 663 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 664 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 665 break; 666 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_LOWER: 667 /* [A00 0; A10 S], suitable for left preconditioning */ 668 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 669 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 670 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 671 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 672 ierr = VecScale(ilinkD->x,-1.);CHKERRQ(ierr); 673 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 674 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 675 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 676 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 677 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 678 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 679 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 680 break; 681 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_UPPER: 682 /* [A00 A01; 0 S], suitable for right preconditioning */ 683 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 684 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 685 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 686 ierr = MatMult(jac->B,ilinkD->y,ilinkA->x);CHKERRQ(ierr); 687 ierr = VecScale(ilinkA->x,-1.);CHKERRQ(ierr); 688 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 689 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 690 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 691 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 692 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 693 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 694 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 695 break; 696 case PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL: 697 /* [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 */ 698 ierr = VecScatterBegin(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 699 ierr = VecScatterEnd(ilinkA->sctx,x,ilinkA->x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 700 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 701 ierr = MatMult(jac->C,ilinkA->y,ilinkD->x);CHKERRQ(ierr); 702 ierr = VecScale(ilinkD->x,-1.0);CHKERRQ(ierr); 703 ierr = VecScatterBegin(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 704 ierr = VecScatterEnd(ilinkD->sctx,x,ilinkD->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 705 706 ierr = KSPSolve(jac->kspschur,ilinkD->x,ilinkD->y);CHKERRQ(ierr); 707 ierr = VecScatterBegin(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 708 ierr = VecScatterEnd(ilinkD->sctx,ilinkD->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 709 710 ierr = MatMult(jac->B,ilinkD->y,ilinkA->y);CHKERRQ(ierr); 711 ierr = VecAXPY(ilinkA->x,-1.0,ilinkA->y);CHKERRQ(ierr); 712 ierr = KSPSolve(ksp,ilinkA->x,ilinkA->y);CHKERRQ(ierr); 713 ierr = VecScatterBegin(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 714 ierr = VecScatterEnd(ilinkA->sctx,ilinkA->y,y,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 715 } 716 PetscFunctionReturn(0); 717 } 718 719 #undef __FUNCT__ 720 #define __FUNCT__ "PCApply_FieldSplit" 721 static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y) 722 { 723 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 724 PetscErrorCode ierr; 725 PC_FieldSplitLink ilink = jac->head; 726 PetscInt cnt; 727 728 PetscFunctionBegin; 729 CHKMEMQ; 730 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 731 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 732 733 if (jac->type == PC_COMPOSITE_ADDITIVE) { 734 if (jac->defaultsplit) { 735 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 736 while (ilink) { 737 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 738 ilink = ilink->next; 739 } 740 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 741 } else { 742 ierr = VecSet(y,0.0);CHKERRQ(ierr); 743 while (ilink) { 744 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 745 ilink = ilink->next; 746 } 747 } 748 } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 749 if (!jac->w1) { 750 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 751 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 752 } 753 ierr = VecSet(y,0.0);CHKERRQ(ierr); 754 ierr = FieldSplitSplitSolveAdd(ilink,x,y);CHKERRQ(ierr); 755 cnt = 1; 756 while (ilink->next) { 757 ilink = ilink->next; 758 /* compute the residual only over the part of the vector needed */ 759 ierr = MatMult(jac->Afield[cnt++],y,ilink->x);CHKERRQ(ierr); 760 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 761 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 762 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 763 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 764 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 765 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 766 } 767 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 768 cnt -= 2; 769 while (ilink->previous) { 770 ilink = ilink->previous; 771 /* compute the residual only over the part of the vector needed */ 772 ierr = MatMult(jac->Afield[cnt--],y,ilink->x);CHKERRQ(ierr); 773 ierr = VecScale(ilink->x,-1.0);CHKERRQ(ierr); 774 ierr = VecScatterBegin(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 775 ierr = VecScatterEnd(ilink->sctx,x,ilink->x,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 776 ierr = KSPSolve(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 777 ierr = VecScatterBegin(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 778 ierr = VecScatterEnd(ilink->sctx,ilink->y,y,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 779 } 780 } 781 } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Unsupported or unknown composition",(int) jac->type); 782 CHKMEMQ; 783 PetscFunctionReturn(0); 784 } 785 786 #define FieldSplitSplitSolveAddTranspose(ilink,xx,yy) \ 787 (VecScatterBegin(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 788 VecScatterEnd(ilink->sctx,xx,ilink->y,INSERT_VALUES,SCATTER_FORWARD) || \ 789 KSPSolveTranspose(ilink->ksp,ilink->y,ilink->x) || \ 790 VecScatterBegin(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE) || \ 791 VecScatterEnd(ilink->sctx,ilink->x,yy,ADD_VALUES,SCATTER_REVERSE)) 792 793 #undef __FUNCT__ 794 #define __FUNCT__ "PCApplyTranspose_FieldSplit" 795 static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc,Vec x,Vec y) 796 { 797 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 798 PetscErrorCode ierr; 799 PC_FieldSplitLink ilink = jac->head; 800 801 PetscFunctionBegin; 802 CHKMEMQ; 803 ierr = VecSetBlockSize(x,jac->bs);CHKERRQ(ierr); 804 ierr = VecSetBlockSize(y,jac->bs);CHKERRQ(ierr); 805 806 if (jac->type == PC_COMPOSITE_ADDITIVE) { 807 if (jac->defaultsplit) { 808 ierr = VecStrideGatherAll(x,jac->x,INSERT_VALUES);CHKERRQ(ierr); 809 while (ilink) { 810 ierr = KSPSolveTranspose(ilink->ksp,ilink->x,ilink->y);CHKERRQ(ierr); 811 ilink = ilink->next; 812 } 813 ierr = VecStrideScatterAll(jac->y,y,INSERT_VALUES);CHKERRQ(ierr); 814 } else { 815 ierr = VecSet(y,0.0);CHKERRQ(ierr); 816 while (ilink) { 817 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 818 ilink = ilink->next; 819 } 820 } 821 } else { 822 if (!jac->w1) { 823 ierr = VecDuplicate(x,&jac->w1);CHKERRQ(ierr); 824 ierr = VecDuplicate(x,&jac->w2);CHKERRQ(ierr); 825 } 826 ierr = VecSet(y,0.0);CHKERRQ(ierr); 827 if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) { 828 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 829 while (ilink->next) { 830 ilink = ilink->next; 831 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 832 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 833 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 834 } 835 while (ilink->previous) { 836 ilink = ilink->previous; 837 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 838 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 839 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 840 } 841 } else { 842 while (ilink->next) { /* get to last entry in linked list */ 843 ilink = ilink->next; 844 } 845 ierr = FieldSplitSplitSolveAddTranspose(ilink,x,y);CHKERRQ(ierr); 846 while (ilink->previous) { 847 ilink = ilink->previous; 848 ierr = MatMultTranspose(pc->mat,y,jac->w1);CHKERRQ(ierr); 849 ierr = VecWAXPY(jac->w2,-1.0,jac->w1,x);CHKERRQ(ierr); 850 ierr = FieldSplitSplitSolveAddTranspose(ilink,jac->w2,y);CHKERRQ(ierr); 851 } 852 } 853 } 854 CHKMEMQ; 855 PetscFunctionReturn(0); 856 } 857 858 #undef __FUNCT__ 859 #define __FUNCT__ "PCReset_FieldSplit" 860 static PetscErrorCode PCReset_FieldSplit(PC pc) 861 { 862 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 863 PetscErrorCode ierr; 864 PC_FieldSplitLink ilink = jac->head,next; 865 866 PetscFunctionBegin; 867 while (ilink) { 868 ierr = KSPReset(ilink->ksp);CHKERRQ(ierr); 869 ierr = VecDestroy(&ilink->x);CHKERRQ(ierr); 870 ierr = VecDestroy(&ilink->y);CHKERRQ(ierr); 871 ierr = VecScatterDestroy(&ilink->sctx);CHKERRQ(ierr); 872 ierr = ISDestroy(&ilink->is);CHKERRQ(ierr); 873 next = ilink->next; 874 ilink = next; 875 } 876 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 877 if (jac->mat && jac->mat != jac->pmat) { 878 ierr = MatDestroyMatrices(jac->nsplits,&jac->mat);CHKERRQ(ierr); 879 } else if (jac->mat) { 880 jac->mat = PETSC_NULL; 881 } 882 if (jac->pmat) {ierr = MatDestroyMatrices(jac->nsplits,&jac->pmat);CHKERRQ(ierr);} 883 if (jac->Afield) {ierr = MatDestroyMatrices(jac->nsplits,&jac->Afield);CHKERRQ(ierr);} 884 ierr = VecDestroy(&jac->w1);CHKERRQ(ierr); 885 ierr = VecDestroy(&jac->w2);CHKERRQ(ierr); 886 ierr = MatDestroy(&jac->schur);CHKERRQ(ierr); 887 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 888 ierr = KSPDestroy(&jac->kspschur);CHKERRQ(ierr); 889 ierr = MatDestroy(&jac->B);CHKERRQ(ierr); 890 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 891 jac->reset = PETSC_TRUE; 892 PetscFunctionReturn(0); 893 } 894 895 #undef __FUNCT__ 896 #define __FUNCT__ "PCDestroy_FieldSplit" 897 static PetscErrorCode PCDestroy_FieldSplit(PC pc) 898 { 899 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 900 PetscErrorCode ierr; 901 PC_FieldSplitLink ilink = jac->head,next; 902 903 PetscFunctionBegin; 904 ierr = PCReset_FieldSplit(pc);CHKERRQ(ierr); 905 while (ilink) { 906 ierr = KSPDestroy(&ilink->ksp);CHKERRQ(ierr); 907 next = ilink->next; 908 ierr = PetscFree(ilink->splitname);CHKERRQ(ierr); 909 ierr = PetscFree(ilink->fields);CHKERRQ(ierr); 910 ierr = PetscFree(ilink);CHKERRQ(ierr); 911 ilink = next; 912 } 913 ierr = PetscFree2(jac->x,jac->y);CHKERRQ(ierr); 914 ierr = PetscFree(pc->data);CHKERRQ(ierr); 915 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","",PETSC_NULL);CHKERRQ(ierr); 916 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","",PETSC_NULL);CHKERRQ(ierr); 917 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","",PETSC_NULL);CHKERRQ(ierr); 918 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","",PETSC_NULL);CHKERRQ(ierr); 919 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","",PETSC_NULL);CHKERRQ(ierr); 920 PetscFunctionReturn(0); 921 } 922 923 #undef __FUNCT__ 924 #define __FUNCT__ "PCSetFromOptions_FieldSplit" 925 static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc) 926 { 927 PetscErrorCode ierr; 928 PetscInt bs; 929 PetscBool flg,stokes = PETSC_FALSE; 930 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 931 PCCompositeType ctype; 932 933 PetscFunctionBegin; 934 ierr = PetscOptionsHead("FieldSplit options");CHKERRQ(ierr); 935 ierr = PetscOptionsBool("-pc_fieldsplit_real_diagonal","Use diagonal blocks of the operator","PCFieldSplitSetRealDiagonal",jac->realdiagonal,&jac->realdiagonal,PETSC_NULL);CHKERRQ(ierr); 936 ierr = PetscOptionsInt("-pc_fieldsplit_block_size","Blocksize that defines number of fields","PCFieldSplitSetBlockSize",jac->bs,&bs,&flg);CHKERRQ(ierr); 937 if (flg) { 938 ierr = PCFieldSplitSetBlockSize(pc,bs);CHKERRQ(ierr); 939 } 940 941 ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_fieldsplit_detect_saddle_point",&stokes,PETSC_NULL);CHKERRQ(ierr); 942 if (stokes) { 943 ierr = PCFieldSplitSetType(pc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 944 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_SELF; 945 } 946 947 ierr = PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&ctype,&flg);CHKERRQ(ierr); 948 if (flg) { 949 ierr = PCFieldSplitSetType(pc,ctype);CHKERRQ(ierr); 950 } 951 952 /* Only setup fields once */ 953 if ((jac->bs > 0) && (jac->nsplits == 0)) { 954 /* only allow user to set fields from command line if bs is already known. 955 otherwise user can set them in PCFieldSplitSetDefaults() */ 956 ierr = PCFieldSplitSetRuntimeSplits_Private(pc);CHKERRQ(ierr); 957 if (jac->splitdefined) {ierr = PetscInfo(pc,"Splits defined using the options database\n");CHKERRQ(ierr);} 958 } 959 if (jac->type == PC_COMPOSITE_SCHUR) { 960 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_factorization_type","Factorization to use","None",PCFieldSplitSchurFactorizationTypes,(PetscEnum)jac->schurfactorization,(PetscEnum*)&jac->schurfactorization,PETSC_NULL);CHKERRQ(ierr); 961 ierr = PetscOptionsEnum("-pc_fieldsplit_schur_precondition","How to build preconditioner for Schur complement","PCFieldSplitSchurPrecondition",PCFieldSplitSchurPreTypes,(PetscEnum)jac->schurpre,(PetscEnum*)&jac->schurpre,PETSC_NULL);CHKERRQ(ierr); 962 } 963 ierr = PetscOptionsTail();CHKERRQ(ierr); 964 PetscFunctionReturn(0); 965 } 966 967 /*------------------------------------------------------------------------------------*/ 968 969 EXTERN_C_BEGIN 970 #undef __FUNCT__ 971 #define __FUNCT__ "PCFieldSplitSetFields_FieldSplit" 972 PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 973 { 974 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 975 PetscErrorCode ierr; 976 PC_FieldSplitLink ilink,next = jac->head; 977 char prefix[128]; 978 PetscInt i; 979 980 PetscFunctionBegin; 981 if (jac->splitdefined) { 982 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 983 PetscFunctionReturn(0); 984 } 985 for (i=0; i<n; i++) { 986 if (fields[i] >= jac->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Field %D requested but only %D exist",fields[i],jac->bs); 987 if (fields[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative field %D requested",fields[i]); 988 } 989 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 990 if (splitname) { 991 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 992 } else { 993 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 994 ierr = PetscSNPrintf(ilink->splitname,2,"%s",jac->nsplits);CHKERRQ(ierr); 995 } 996 ierr = PetscMalloc(n*sizeof(PetscInt),&ilink->fields);CHKERRQ(ierr); 997 ierr = PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));CHKERRQ(ierr); 998 ilink->nfields = n; 999 ilink->next = PETSC_NULL; 1000 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1001 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1002 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1003 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1004 1005 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1006 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1007 1008 if (!next) { 1009 jac->head = ilink; 1010 ilink->previous = PETSC_NULL; 1011 } else { 1012 while (next->next) { 1013 next = next->next; 1014 } 1015 next->next = ilink; 1016 ilink->previous = next; 1017 } 1018 jac->nsplits++; 1019 PetscFunctionReturn(0); 1020 } 1021 EXTERN_C_END 1022 1023 EXTERN_C_BEGIN 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit_Schur" 1026 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc,PetscInt *n,KSP **subksp) 1027 { 1028 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1029 PetscErrorCode ierr; 1030 1031 PetscFunctionBegin; 1032 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1033 ierr = MatSchurComplementGetKSP(jac->schur,*subksp);CHKERRQ(ierr); 1034 (*subksp)[1] = jac->kspschur; 1035 if (n) *n = jac->nsplits; 1036 PetscFunctionReturn(0); 1037 } 1038 EXTERN_C_END 1039 1040 EXTERN_C_BEGIN 1041 #undef __FUNCT__ 1042 #define __FUNCT__ "PCFieldSplitGetSubKSP_FieldSplit" 1043 PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp) 1044 { 1045 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1046 PetscErrorCode ierr; 1047 PetscInt cnt = 0; 1048 PC_FieldSplitLink ilink = jac->head; 1049 1050 PetscFunctionBegin; 1051 ierr = PetscMalloc(jac->nsplits*sizeof(KSP),subksp);CHKERRQ(ierr); 1052 while (ilink) { 1053 (*subksp)[cnt++] = ilink->ksp; 1054 ilink = ilink->next; 1055 } 1056 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); 1057 if (n) *n = jac->nsplits; 1058 PetscFunctionReturn(0); 1059 } 1060 EXTERN_C_END 1061 1062 EXTERN_C_BEGIN 1063 #undef __FUNCT__ 1064 #define __FUNCT__ "PCFieldSplitSetIS_FieldSplit" 1065 PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc,const char splitname[],IS is) 1066 { 1067 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1068 PetscErrorCode ierr; 1069 PC_FieldSplitLink ilink, next = jac->head; 1070 char prefix[128]; 1071 1072 PetscFunctionBegin; 1073 if (jac->splitdefined) { 1074 ierr = PetscInfo1(pc,"Ignoring new split \"%s\" because the splits have already been defined\n",splitname);CHKERRQ(ierr); 1075 PetscFunctionReturn(0); 1076 } 1077 ierr = PetscNew(struct _PC_FieldSplitLink,&ilink);CHKERRQ(ierr); 1078 if (splitname) { 1079 ierr = PetscStrallocpy(splitname,&ilink->splitname);CHKERRQ(ierr); 1080 } else { 1081 ierr = PetscMalloc(3*sizeof(char),&ilink->splitname);CHKERRQ(ierr); 1082 ierr = PetscSNPrintf(ilink->splitname,3,"%D",jac->nsplits);CHKERRQ(ierr); 1083 } 1084 ilink->is = is; 1085 ierr = PetscObjectReference((PetscObject)is);CHKERRQ(ierr); 1086 ilink->next = PETSC_NULL; 1087 ierr = KSPCreate(((PetscObject)pc)->comm,&ilink->ksp);CHKERRQ(ierr); 1088 ierr = PetscObjectIncrementTabLevel((PetscObject)ilink->ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1089 ierr = KSPSetType(ilink->ksp,KSPPREONLY);CHKERRQ(ierr); 1090 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)ilink->ksp);CHKERRQ(ierr); 1091 1092 ierr = PetscSNPrintf(prefix,sizeof prefix,"%sfieldsplit_%s_",((PetscObject)pc)->prefix?((PetscObject)pc)->prefix:"",ilink->splitname);CHKERRQ(ierr); 1093 ierr = KSPSetOptionsPrefix(ilink->ksp,prefix);CHKERRQ(ierr); 1094 1095 if (!next) { 1096 jac->head = ilink; 1097 ilink->previous = PETSC_NULL; 1098 } else { 1099 while (next->next) { 1100 next = next->next; 1101 } 1102 next->next = ilink; 1103 ilink->previous = next; 1104 } 1105 jac->nsplits++; 1106 1107 PetscFunctionReturn(0); 1108 } 1109 EXTERN_C_END 1110 1111 #undef __FUNCT__ 1112 #define __FUNCT__ "PCFieldSplitSetFields" 1113 /*@ 1114 PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner 1115 1116 Logically Collective on PC 1117 1118 Input Parameters: 1119 + pc - the preconditioner context 1120 . splitname - name of this split, if PETSC_NULL the number of the split is used 1121 . n - the number of fields in this split 1122 - fields - the fields in this split 1123 1124 Level: intermediate 1125 1126 Notes: Use PCFieldSplitSetIS() to set a completely general set of indices as a field. 1127 1128 The PCFieldSplitSetFields() is for defining fields as a strided blocks. For example, if the block 1129 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 1130 0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x1x3x5x7.. x12x45x78x.... 1131 where the numbered entries indicate what is in the field. 1132 1133 This function is called once per split (it creates a new split each time). Solve options 1134 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1135 1136 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize(), PCFieldSplitSetIS() 1137 1138 @*/ 1139 PetscErrorCode PCFieldSplitSetFields(PC pc,const char splitname[],PetscInt n,const PetscInt *fields) 1140 { 1141 PetscErrorCode ierr; 1142 1143 PetscFunctionBegin; 1144 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1145 PetscValidCharPointer(splitname,2); 1146 if (n < 1) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Provided number of fields %D in split \"%s\" not positive",n,splitname); 1147 PetscValidIntPointer(fields,3); 1148 ierr = PetscTryMethod(pc,"PCFieldSplitSetFields_C",(PC,const char[],PetscInt,const PetscInt *),(pc,splitname,n,fields));CHKERRQ(ierr); 1149 PetscFunctionReturn(0); 1150 } 1151 1152 #undef __FUNCT__ 1153 #define __FUNCT__ "PCFieldSplitSetIS" 1154 /*@ 1155 PCFieldSplitSetIS - Sets the exact elements for field 1156 1157 Logically Collective on PC 1158 1159 Input Parameters: 1160 + pc - the preconditioner context 1161 . splitname - name of this split, if PETSC_NULL the number of the split is used 1162 - is - the index set that defines the vector elements in this field 1163 1164 1165 Notes: 1166 Use PCFieldSplitSetFields(), for fields defined by strided types. 1167 1168 This function is called once per split (it creates a new split each time). Solve options 1169 for this split will be available under the prefix -fieldsplit_SPLITNAME_. 1170 1171 Level: intermediate 1172 1173 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetBlockSize() 1174 1175 @*/ 1176 PetscErrorCode PCFieldSplitSetIS(PC pc,const char splitname[],IS is) 1177 { 1178 PetscErrorCode ierr; 1179 1180 PetscFunctionBegin; 1181 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1182 if (splitname) PetscValidCharPointer(splitname,2); 1183 PetscValidHeaderSpecific(is,IS_CLASSID,3); 1184 ierr = PetscTryMethod(pc,"PCFieldSplitSetIS_C",(PC,const char[],IS),(pc,splitname,is));CHKERRQ(ierr); 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "PCFieldSplitGetIS" 1190 /*@ 1191 PCFieldSplitGetIS - Retrieves the elements for a field as an IS 1192 1193 Logically Collective on PC 1194 1195 Input Parameters: 1196 + pc - the preconditioner context 1197 - splitname - name of this split 1198 1199 Output Parameter: 1200 - is - the index set that defines the vector elements in this field, or PETSC_NULL if the field is not found 1201 1202 Level: intermediate 1203 1204 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetIS() 1205 1206 @*/ 1207 PetscErrorCode PCFieldSplitGetIS(PC pc,const char splitname[],IS *is) 1208 { 1209 PetscErrorCode ierr; 1210 1211 PetscFunctionBegin; 1212 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1213 PetscValidCharPointer(splitname,2); 1214 PetscValidPointer(is,3); 1215 { 1216 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1217 PC_FieldSplitLink ilink = jac->head; 1218 PetscBool found; 1219 1220 *is = PETSC_NULL; 1221 while(ilink) { 1222 ierr = PetscStrcmp(ilink->splitname, splitname, &found);CHKERRQ(ierr); 1223 if (found) { 1224 *is = ilink->is; 1225 break; 1226 } 1227 ilink = ilink->next; 1228 } 1229 } 1230 PetscFunctionReturn(0); 1231 } 1232 1233 #undef __FUNCT__ 1234 #define __FUNCT__ "PCFieldSplitSetBlockSize" 1235 /*@ 1236 PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the 1237 fieldsplit preconditioner. If not set the matrix block size is used. 1238 1239 Logically Collective on PC 1240 1241 Input Parameters: 1242 + pc - the preconditioner context 1243 - bs - the block size 1244 1245 Level: intermediate 1246 1247 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields() 1248 1249 @*/ 1250 PetscErrorCode PCFieldSplitSetBlockSize(PC pc,PetscInt bs) 1251 { 1252 PetscErrorCode ierr; 1253 1254 PetscFunctionBegin; 1255 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1256 PetscValidLogicalCollectiveInt(pc,bs,2); 1257 ierr = PetscTryMethod(pc,"PCFieldSplitSetBlockSize_C",(PC,PetscInt),(pc,bs));CHKERRQ(ierr); 1258 PetscFunctionReturn(0); 1259 } 1260 1261 #undef __FUNCT__ 1262 #define __FUNCT__ "PCFieldSplitGetSubKSP" 1263 /*@C 1264 PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits 1265 1266 Collective on KSP 1267 1268 Input Parameter: 1269 . pc - the preconditioner context 1270 1271 Output Parameters: 1272 + n - the number of splits 1273 - pc - the array of KSP contexts 1274 1275 Note: 1276 After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed by the user 1277 (not the KSP just the array that contains them). 1278 1279 You must call KSPSetUp() before calling PCFieldSplitGetSubKSP(). 1280 1281 Level: advanced 1282 1283 .seealso: PCFIELDSPLIT 1284 @*/ 1285 PetscErrorCode PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[]) 1286 { 1287 PetscErrorCode ierr; 1288 1289 PetscFunctionBegin; 1290 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1291 if (n) PetscValidIntPointer(n,2); 1292 ierr = PetscUseMethod(pc,"PCFieldSplitGetSubKSP_C",(PC,PetscInt*,KSP **),(pc,n,subksp));CHKERRQ(ierr); 1293 PetscFunctionReturn(0); 1294 } 1295 1296 #undef __FUNCT__ 1297 #define __FUNCT__ "PCFieldSplitSchurPrecondition" 1298 /*@ 1299 PCFieldSplitSchurPrecondition - Indicates if the Schur complement is preconditioned by a preconditioner constructed by the 1300 A11 matrix. Otherwise no preconditioner is used. 1301 1302 Collective on PC 1303 1304 Input Parameters: 1305 + pc - the preconditioner context 1306 . ptype - which matrix to use for preconditioning the Schur complement, PC_FIELDSPLIT_SCHUR_PRE_DIAG (diag) is default 1307 - userpre - matrix to use for preconditioning, or PETSC_NULL 1308 1309 Options Database: 1310 . -pc_fieldsplit_schur_precondition <self,user,diag> default is diag 1311 1312 Notes: 1313 $ If ptype is 1314 $ user then the preconditioner for the Schur complement is generated by the provided matrix (pre argument 1315 $ to this function). 1316 $ diag then the preconditioner for the Schur complement is generated by the block diagonal part of the original 1317 $ matrix associated with the Schur complement (i.e. A11) 1318 $ self the preconditioner for the Schur complement is generated from the Schur complement matrix itself: 1319 $ The only preconditioner that currently works directly with the Schur complement matrix object is the PCLSC 1320 $ preconditioner 1321 1322 When solving a saddle point problem, where the A11 block is identically zero, using diag as the ptype only makes sense 1323 with the additional option -fieldsplit_1_pc_type none. Usually for saddle point problems one would use a ptype of self and 1324 -fieldsplit_1_pc_type lsc which uses the least squares commutator compute a preconditioner for the Schur complement. 1325 1326 Developer Notes: This is a terrible name, gives no good indication of what the function does and should also have Set in 1327 the name since it sets a proceedure to use. 1328 1329 Level: intermediate 1330 1331 .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT, PCFieldSplitSetFields(), PCFieldSplitSchurPreType, PCLSC 1332 1333 @*/ 1334 PetscErrorCode PCFieldSplitSchurPrecondition(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1335 { 1336 PetscErrorCode ierr; 1337 1338 PetscFunctionBegin; 1339 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1340 ierr = PetscTryMethod(pc,"PCFieldSplitSchurPrecondition_C",(PC,PCFieldSplitSchurPreType,Mat),(pc,ptype,pre));CHKERRQ(ierr); 1341 PetscFunctionReturn(0); 1342 } 1343 1344 EXTERN_C_BEGIN 1345 #undef __FUNCT__ 1346 #define __FUNCT__ "PCFieldSplitSchurPrecondition_FieldSplit" 1347 PetscErrorCode PCFieldSplitSchurPrecondition_FieldSplit(PC pc,PCFieldSplitSchurPreType ptype,Mat pre) 1348 { 1349 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1350 PetscErrorCode ierr; 1351 1352 PetscFunctionBegin; 1353 jac->schurpre = ptype; 1354 if (pre) { 1355 ierr = MatDestroy(&jac->schur_user);CHKERRQ(ierr); 1356 jac->schur_user = pre; 1357 ierr = PetscObjectReference((PetscObject)jac->schur_user);CHKERRQ(ierr); 1358 } 1359 PetscFunctionReturn(0); 1360 } 1361 EXTERN_C_END 1362 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "PCFieldSplitGetSchurBlocks" 1365 /*@C 1366 PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement 1367 1368 Collective on KSP 1369 1370 Input Parameter: 1371 . pc - the preconditioner context 1372 1373 Output Parameters: 1374 + A00 - the (0,0) block 1375 . A01 - the (0,1) block 1376 . A10 - the (1,0) block 1377 - A11 - the (1,1) block 1378 1379 Level: advanced 1380 1381 .seealso: PCFIELDSPLIT 1382 @*/ 1383 PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc,Mat *A00,Mat *A01,Mat *A10, Mat *A11) 1384 { 1385 PC_FieldSplit *jac = (PC_FieldSplit *) pc->data; 1386 1387 PetscFunctionBegin; 1388 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1389 if (jac->type != PC_COMPOSITE_SCHUR) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach."); 1390 if (A00) *A00 = jac->pmat[0]; 1391 if (A01) *A01 = jac->B; 1392 if (A10) *A10 = jac->C; 1393 if (A11) *A11 = jac->pmat[1]; 1394 PetscFunctionReturn(0); 1395 } 1396 1397 EXTERN_C_BEGIN 1398 #undef __FUNCT__ 1399 #define __FUNCT__ "PCFieldSplitSetType_FieldSplit" 1400 PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type) 1401 { 1402 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1403 PetscErrorCode ierr; 1404 1405 PetscFunctionBegin; 1406 jac->type = type; 1407 if (type == PC_COMPOSITE_SCHUR) { 1408 pc->ops->apply = PCApply_FieldSplit_Schur; 1409 pc->ops->view = PCView_FieldSplit_Schur; 1410 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit_Schur",PCFieldSplitGetSubKSP_FieldSplit_Schur);CHKERRQ(ierr); 1411 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","PCFieldSplitSchurPrecondition_FieldSplit",PCFieldSplitSchurPrecondition_FieldSplit);CHKERRQ(ierr); 1412 1413 } else { 1414 pc->ops->apply = PCApply_FieldSplit; 1415 pc->ops->view = PCView_FieldSplit; 1416 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1417 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSchurPrecondition_C","",0);CHKERRQ(ierr); 1418 } 1419 PetscFunctionReturn(0); 1420 } 1421 EXTERN_C_END 1422 1423 EXTERN_C_BEGIN 1424 #undef __FUNCT__ 1425 #define __FUNCT__ "PCFieldSplitSetBlockSize_FieldSplit" 1426 PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc,PetscInt bs) 1427 { 1428 PC_FieldSplit *jac = (PC_FieldSplit*)pc->data; 1429 1430 PetscFunctionBegin; 1431 if (bs < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Blocksize must be positive, you gave %D",bs); 1432 if (jac->bs > 0 && jac->bs != bs) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot change fieldsplit blocksize from %D to %D after it has been set",jac->bs,bs); 1433 jac->bs = bs; 1434 PetscFunctionReturn(0); 1435 } 1436 EXTERN_C_END 1437 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "PCFieldSplitSetType" 1440 /*@ 1441 PCFieldSplitSetType - Sets the type of fieldsplit preconditioner. 1442 1443 Collective on PC 1444 1445 Input Parameter: 1446 . pc - the preconditioner context 1447 . type - PC_COMPOSITE_ADDITIVE, PC_COMPOSITE_MULTIPLICATIVE (default), PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL, PC_COMPOSITE_SCHUR 1448 1449 Options Database Key: 1450 . -pc_fieldsplit_type <type: one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type 1451 1452 Level: Developer 1453 1454 .keywords: PC, set, type, composite preconditioner, additive, multiplicative 1455 1456 .seealso: PCCompositeSetType() 1457 1458 @*/ 1459 PetscErrorCode PCFieldSplitSetType(PC pc,PCCompositeType type) 1460 { 1461 PetscErrorCode ierr; 1462 1463 PetscFunctionBegin; 1464 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1465 ierr = PetscTryMethod(pc,"PCFieldSplitSetType_C",(PC,PCCompositeType),(pc,type));CHKERRQ(ierr); 1466 PetscFunctionReturn(0); 1467 } 1468 1469 /* -------------------------------------------------------------------------------------*/ 1470 /*MC 1471 PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual 1472 fields or groups of fields. See the users manual section "Solving Block Matrices" for more details. 1473 1474 To set options on the solvers for each block append -fieldsplit_ to all the PC 1475 options database keys. For example, -fieldsplit_pc_type ilu -fieldsplit_pc_factor_levels 1 1476 1477 To set the options on the solvers separate for each block call PCFieldSplitGetSubKSP() 1478 and set the options directly on the resulting KSP object 1479 1480 Level: intermediate 1481 1482 Options Database Keys: 1483 + -pc_fieldsplit_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split 1484 . -pc_fieldsplit_default - automatically add any fields to additional splits that have not 1485 been supplied explicitly by -pc_fieldsplit_%d_fields 1486 . -pc_fieldsplit_block_size <bs> - size of block that defines fields (i.e. there are bs fields) 1487 . -pc_fieldsplit_type <additive,multiplicative,schur,symmetric_multiplicative> 1488 . -pc_fieldsplit_schur_precondition <true,false> default is true 1489 . -pc_fieldsplit_detect_saddle_point - automatically finds rows with zero or negative diagonal and uses Schur complement with no preconditioner as the solver 1490 . -fieldsplit_NAME_ksp_* - control inner linear solver, NAME is a sequential integer if unspecified, otherwise use name provided in PCFieldSplitSetIS() or the name associated with the field in the DM passed to PCSetDM() (or a higher level function) 1491 - -fieldsplit_NAME_pc_* - control inner preconditioner, NAME has same semantics as above 1492 1493 Notes: 1494 Use PCFieldSplitSetFields() to set fields defined by "strided" entries and PCFieldSplitSetIS() 1495 to define a field by an arbitrary collection of entries. 1496 1497 If no fields are set the default is used. The fields are defined by entries strided by bs, 1498 beginning at 0 then 1, etc to bs-1. The block size can be set with PCFieldSplitSetBlockSize(), 1499 if this is not called the block size defaults to the blocksize of the second matrix passed 1500 to KSPSetOperators()/PCSetOperators(). 1501 1502 For the Schur complement preconditioner if J = ( A00 A01 ) 1503 ( A10 A11 ) 1504 the preconditioner is 1505 (I -B ksp(A00)) ( inv(A00) 0 (I 0 ) 1506 (0 I ) ( 0 ksp(S) ) (-A10 ksp(A00) I ) 1507 where the action of inv(A00) is applied using the KSP solver with prefix -fieldsplit_0_. The action of 1508 ksp(S) is computed using the KSP solver with prefix -fieldsplit_splitname_ (where splitname was given in providing the SECOND split or 1 if not give). 1509 For PCFieldSplitGetKSP() when field number is 1510 0 it returns the KSP associated with -fieldsplit_0_ while field number 1 gives -fieldsplit_1_ KSP. By default 1511 A11 is used to construct a preconditioner for S, use PCFieldSplitSchurPrecondition() to turn on or off this 1512 option. You can use the preconditioner PCLSC to precondition the Schur complement with -fieldsplit_1_pc_type lsc 1513 1514 If only one set of indices (one IS) is provided with PCFieldSplitSetIS() then the complement of that IS 1515 is used automatically for a second block. 1516 1517 The fieldsplit preconditioner cannot currently be used with the BAIJ or SBAIJ data formats if the blocksize is larger than 1. 1518 Generally it should be used with the AIJ format. 1519 1520 The forms of these preconditioners are closely related if not identical to forms derived as "Distributive Iterations", see, 1521 for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling. Note that one can also use PCFIELDSPLIT 1522 inside a smoother resulting in "Distributive Smoothers". 1523 1524 Concepts: physics based preconditioners, block preconditioners 1525 1526 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCLSC, 1527 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition() 1528 M*/ 1529 1530 EXTERN_C_BEGIN 1531 #undef __FUNCT__ 1532 #define __FUNCT__ "PCCreate_FieldSplit" 1533 PetscErrorCode PCCreate_FieldSplit(PC pc) 1534 { 1535 PetscErrorCode ierr; 1536 PC_FieldSplit *jac; 1537 1538 PetscFunctionBegin; 1539 ierr = PetscNewLog(pc,PC_FieldSplit,&jac);CHKERRQ(ierr); 1540 jac->bs = -1; 1541 jac->nsplits = 0; 1542 jac->type = PC_COMPOSITE_MULTIPLICATIVE; 1543 jac->schurpre = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */ 1544 jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACTORIZATION_FULL; 1545 1546 pc->data = (void*)jac; 1547 1548 pc->ops->apply = PCApply_FieldSplit; 1549 pc->ops->applytranspose = PCApplyTranspose_FieldSplit; 1550 pc->ops->setup = PCSetUp_FieldSplit; 1551 pc->ops->reset = PCReset_FieldSplit; 1552 pc->ops->destroy = PCDestroy_FieldSplit; 1553 pc->ops->setfromoptions = PCSetFromOptions_FieldSplit; 1554 pc->ops->view = PCView_FieldSplit; 1555 pc->ops->applyrichardson = 0; 1556 1557 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit", 1558 PCFieldSplitGetSubKSP_FieldSplit);CHKERRQ(ierr); 1559 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit", 1560 PCFieldSplitSetFields_FieldSplit);CHKERRQ(ierr); 1561 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetIS_C","PCFieldSplitSetIS_FieldSplit", 1562 PCFieldSplitSetIS_FieldSplit);CHKERRQ(ierr); 1563 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit", 1564 PCFieldSplitSetType_FieldSplit);CHKERRQ(ierr); 1565 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetBlockSize_C","PCFieldSplitSetBlockSize_FieldSplit", 1566 PCFieldSplitSetBlockSize_FieldSplit);CHKERRQ(ierr); 1567 PetscFunctionReturn(0); 1568 } 1569 EXTERN_C_END 1570 1571 1572 1573