1 2 #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/ 3 4 typedef struct { 5 PetscBool allocated; 6 PetscBool scalediag; 7 KSP kspL; 8 Vec scale; 9 Vec x0,y0,x1; 10 Mat L; /* keep a copy to reuse when obtained with L = A10*A01 */ 11 } PC_LSC; 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "PCLSCAllocate_Private" 15 static PetscErrorCode PCLSCAllocate_Private(PC pc) 16 { 17 PC_LSC *lsc = (PC_LSC*)pc->data; 18 Mat A; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 if (lsc->allocated) PetscFunctionReturn(0); 23 ierr = KSPCreate(PetscObjectComm((PetscObject)pc),&lsc->kspL);CHKERRQ(ierr); 24 ierr = PetscObjectIncrementTabLevel((PetscObject)lsc->kspL,(PetscObject)pc,1);CHKERRQ(ierr); 25 ierr = KSPSetType(lsc->kspL,KSPPREONLY);CHKERRQ(ierr); 26 ierr = KSPSetOptionsPrefix(lsc->kspL,((PetscObject)pc)->prefix);CHKERRQ(ierr); 27 ierr = KSPAppendOptionsPrefix(lsc->kspL,"lsc_");CHKERRQ(ierr); 28 ierr = KSPSetFromOptions(lsc->kspL);CHKERRQ(ierr); 29 ierr = MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,NULL,NULL,NULL);CHKERRQ(ierr); 30 ierr = MatCreateVecs(A,&lsc->x0,&lsc->y0);CHKERRQ(ierr); 31 ierr = MatCreateVecs(pc->pmat,&lsc->x1,NULL);CHKERRQ(ierr); 32 if (lsc->scalediag) { 33 ierr = VecDuplicate(lsc->x0,&lsc->scale);CHKERRQ(ierr); 34 } 35 lsc->allocated = PETSC_TRUE; 36 PetscFunctionReturn(0); 37 } 38 39 #undef __FUNCT__ 40 #define __FUNCT__ "PCSetUp_LSC" 41 static PetscErrorCode PCSetUp_LSC(PC pc) 42 { 43 PC_LSC *lsc = (PC_LSC*)pc->data; 44 Mat L,Lp,B,C; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 ierr = PCLSCAllocate_Private(pc);CHKERRQ(ierr); 49 ierr = PetscObjectQuery((PetscObject)pc->mat,"LSC_L",(PetscObject*)&L);CHKERRQ(ierr); 50 if (!L) {ierr = PetscObjectQuery((PetscObject)pc->pmat,"LSC_L",(PetscObject*)&L);CHKERRQ(ierr);} 51 ierr = PetscObjectQuery((PetscObject)pc->pmat,"LSC_Lp",(PetscObject*)&Lp);CHKERRQ(ierr); 52 if (!Lp) {ierr = PetscObjectQuery((PetscObject)pc->mat,"LSC_Lp",(PetscObject*)&Lp);CHKERRQ(ierr);} 53 if (!L) { 54 ierr = MatSchurComplementGetSubMatrices(pc->mat,NULL,NULL,&B,&C,NULL);CHKERRQ(ierr); 55 if (!lsc->L) { 56 ierr = MatMatMult(C,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&lsc->L);CHKERRQ(ierr); 57 } else { 58 ierr = MatMatMult(C,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&lsc->L);CHKERRQ(ierr); 59 } 60 Lp = L = lsc->L; 61 } 62 if (lsc->scale) { 63 Mat Ap; 64 ierr = MatSchurComplementGetSubMatrices(pc->mat,NULL,&Ap,NULL,NULL,NULL);CHKERRQ(ierr); 65 ierr = MatGetDiagonal(Ap,lsc->scale);CHKERRQ(ierr); /* Should be the mass matrix, but we don't have plumbing for that yet */ 66 ierr = VecReciprocal(lsc->scale);CHKERRQ(ierr); 67 } 68 ierr = KSPSetOperators(lsc->kspL,L,Lp);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "PCApply_LSC" 74 static PetscErrorCode PCApply_LSC(PC pc,Vec x,Vec y) 75 { 76 PC_LSC *lsc = (PC_LSC*)pc->data; 77 Mat A,B,C; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = MatSchurComplementGetSubMatrices(pc->mat,&A,NULL,&B,&C,NULL);CHKERRQ(ierr); 82 ierr = KSPSolve(lsc->kspL,x,lsc->x1);CHKERRQ(ierr); 83 ierr = MatMult(B,lsc->x1,lsc->x0);CHKERRQ(ierr); 84 if (lsc->scale) { 85 ierr = VecPointwiseMult(lsc->x0,lsc->x0,lsc->scale);CHKERRQ(ierr); 86 } 87 ierr = MatMult(A,lsc->x0,lsc->y0);CHKERRQ(ierr); 88 if (lsc->scale) { 89 ierr = VecPointwiseMult(lsc->y0,lsc->y0,lsc->scale);CHKERRQ(ierr); 90 } 91 ierr = MatMult(C,lsc->y0,lsc->x1);CHKERRQ(ierr); 92 ierr = KSPSolve(lsc->kspL,lsc->x1,y);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "PCReset_LSC" 98 static PetscErrorCode PCReset_LSC(PC pc) 99 { 100 PC_LSC *lsc = (PC_LSC*)pc->data; 101 PetscErrorCode ierr; 102 103 PetscFunctionBegin; 104 ierr = VecDestroy(&lsc->x0);CHKERRQ(ierr); 105 ierr = VecDestroy(&lsc->y0);CHKERRQ(ierr); 106 ierr = VecDestroy(&lsc->x1);CHKERRQ(ierr); 107 ierr = VecDestroy(&lsc->scale);CHKERRQ(ierr); 108 ierr = KSPDestroy(&lsc->kspL);CHKERRQ(ierr); 109 ierr = MatDestroy(&lsc->L);CHKERRQ(ierr); 110 PetscFunctionReturn(0); 111 } 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "PCDestroy_LSC" 115 static PetscErrorCode PCDestroy_LSC(PC pc) 116 { 117 PetscErrorCode ierr; 118 119 PetscFunctionBegin; 120 ierr = PCReset_LSC(pc);CHKERRQ(ierr); 121 ierr = PetscFree(pc->data);CHKERRQ(ierr); 122 PetscFunctionReturn(0); 123 } 124 125 #undef __FUNCT__ 126 #define __FUNCT__ "PCSetFromOptions_LSC" 127 static PetscErrorCode PCSetFromOptions_LSC(PetscOptionsObjectType *PetscOptionsObject,PC pc) 128 { 129 PC_LSC *lsc = (PC_LSC*)pc->data; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = PetscOptionsHead(PetscOptionsObject,"LSC options");CHKERRQ(ierr); 134 { 135 ierr = PetscOptionsBool("-pc_lsc_scale_diag","Use diagonal of velocity block (A) for scaling","None",lsc->scalediag,&lsc->scalediag,NULL);CHKERRQ(ierr); 136 ierr = PetscOptionsTail();CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 #undef __FUNCT__ 141 #define __FUNCT__ "PCView_LSC" 142 static PetscErrorCode PCView_LSC(PC pc,PetscViewer viewer) 143 { 144 PC_LSC *jac = (PC_LSC*)pc->data; 145 PetscErrorCode ierr; 146 PetscBool iascii; 147 148 PetscFunctionBegin; 149 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 150 if (iascii) { 151 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 152 ierr = KSPView(jac->kspL,viewer);CHKERRQ(ierr); 153 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 154 } 155 PetscFunctionReturn(0); 156 } 157 158 /*MC 159 PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators 160 161 Options Database Key: 162 . -pc_lsc_scale_diag - Use the diagonal of A for scaling 163 164 Level: intermediate 165 166 Notes: 167 This preconditioner will normally be used with PCFieldSplit to precondition the Schur complement, but 168 it can be used for any Schur complement system. Consider the Schur complement 169 170 .vb 171 S = A11 - A10 inv(A00) A01 172 .ve 173 174 PCLSC currently doesn't do anything with A11, so let's assume it is 0. The idea is that a good approximation to 175 inv(S) is given by 176 177 .vb 178 inv(A10 A01) A10 A00 A01 inv(A10 A01) 179 .ve 180 181 The product A10 A01 can be computed for you, but you can provide it (this is 182 usually more efficient anyway). In the case of incompressible flow, A10 A10 is a Laplacian, call it L. The current 183 interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix. 184 185 If you had called KSPSetOperators(ksp,S,Sp), S should have type MATSCHURCOMPLEMENT and Sp can be any type you 186 like (PCLSC doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp". 187 For example, you might have setup code like this 188 189 .vb 190 PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L); 191 PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp); 192 .ve 193 194 And then your Jacobian assembly would look like 195 196 .vb 197 PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L); 198 PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp); 199 if (L) { assembly L } 200 if (Lp) { assemble Lp } 201 .ve 202 203 With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L 204 205 .vb 206 -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml 207 .ve 208 209 Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners. 210 211 References: 212 + Elman, Howle, Shadid, Shuttleworth, and Tuminaro, Block preconditioners based on approximate commutators, 2006. 213 - Silvester, Elman, Kay, Wathen, Efficient preconditioning of the linearized Navier-Stokes equations for incompressible flow, 2001. 214 215 Concepts: physics based preconditioners, block preconditioners 216 217 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCFIELDSPLIT, 218 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSetSchurPre(), 219 MatCreateSchurComplement() 220 M*/ 221 222 #undef __FUNCT__ 223 #define __FUNCT__ "PCCreate_LSC" 224 PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC pc) 225 { 226 PC_LSC *lsc; 227 PetscErrorCode ierr; 228 229 PetscFunctionBegin; 230 ierr = PetscNewLog(pc,&lsc);CHKERRQ(ierr); 231 pc->data = (void*)lsc; 232 233 pc->ops->apply = PCApply_LSC; 234 pc->ops->applytranspose = 0; 235 pc->ops->setup = PCSetUp_LSC; 236 pc->ops->reset = PCReset_LSC; 237 pc->ops->destroy = PCDestroy_LSC; 238 pc->ops->setfromoptions = PCSetFromOptions_LSC; 239 pc->ops->view = PCView_LSC; 240 pc->ops->applyrichardson = 0; 241 PetscFunctionReturn(0); 242 } 243