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