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