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