1 #define PETSCKSP_DLL 2 3 4 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 5 6 typedef struct { 7 PetscTruth allocated; 8 PetscTruth scalediag; 9 KSP kspL; 10 Vec scale; 11 Vec x0,y0,x1; 12 } PC_LSC; 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "PCLSCAllocate_Private" 16 static PetscErrorCode PCLSCAllocate_Private(PC pc) 17 { 18 PC_LSC *lsc = (PC_LSC*)pc->data; 19 Mat A; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 if (lsc->allocated) PetscFunctionReturn(0); 24 ierr = KSPCreate(((PetscObject)pc)->comm,&lsc->kspL);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; 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 (lsc->scale) { 52 Mat Ap; 53 ierr = MatSchurComplementGetSubmatrices(pc->mat,PETSC_NULL,&Ap,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 54 ierr = MatGetDiagonal(Ap,lsc->scale);CHKERRQ(ierr); /* Should be the mass matrix, but we don't have plumbing for that yet */ 55 ierr = VecReciprocal(lsc->scale);CHKERRQ(ierr); 56 } 57 ierr = KSPSetOperators(lsc->kspL,L,Lp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "PCApply_LSC" 63 static PetscErrorCode PCApply_LSC(PC pc,Vec x,Vec y) 64 { 65 PC_LSC *lsc = (PC_LSC*)pc->data; 66 Mat A,B,C; 67 PetscErrorCode ierr; 68 69 PetscFunctionBegin; 70 ierr = MatSchurComplementGetSubmatrices(pc->mat,&A,PETSC_NULL,&B,&C,PETSC_NULL);CHKERRQ(ierr); 71 ierr = KSPSolve(lsc->kspL,x,lsc->x1);CHKERRQ(ierr); 72 ierr = MatMult(B,lsc->x1,lsc->x0);CHKERRQ(ierr); 73 if (lsc->scale) { 74 ierr = VecPointwiseMult(lsc->x0,lsc->x0,lsc->scale);CHKERRQ(ierr); 75 } 76 ierr = MatMult(A,lsc->x0,lsc->y0);CHKERRQ(ierr); 77 if (lsc->scale) { 78 ierr = VecPointwiseMult(lsc->y0,lsc->y0,lsc->scale);CHKERRQ(ierr); 79 } 80 ierr = MatMult(C,lsc->y0,lsc->x1);CHKERRQ(ierr); 81 ierr = KSPSolve(lsc->kspL,lsc->x1,y);CHKERRQ(ierr); 82 PetscFunctionReturn(0); 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "PCDestroy_LSC" 87 static PetscErrorCode PCDestroy_LSC(PC pc) 88 { 89 PC_LSC *lsc = (PC_LSC*)pc->data; 90 PetscErrorCode ierr; 91 92 PetscFunctionBegin; 93 if (lsc->x0) {ierr = VecDestroy(lsc->x0);CHKERRQ(ierr);} 94 if (lsc->y0) {ierr = VecDestroy(lsc->y0);CHKERRQ(ierr);} 95 if (lsc->x1) {ierr = VecDestroy(lsc->x1);CHKERRQ(ierr);} 96 if (lsc->scale) {ierr = VecDestroy(lsc->scale);CHKERRQ(ierr);} 97 if (lsc->kspL) {ierr = KSPDestroy(lsc->kspL);CHKERRQ(ierr);} 98 ierr = PetscFree(lsc);CHKERRQ(ierr); 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "PCSetFromOptions_LSC" 104 static PetscErrorCode PCSetFromOptions_LSC(PC pc) 105 { 106 PC_LSC *lsc = (PC_LSC*)pc->data; 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 ierr = PetscOptionsHead("LSC options");CHKERRQ(ierr); 111 { 112 ierr = PetscOptionsTruth("-pc_lsc_scale_diag","Use diagonal of velocity block (A) for scaling","None",lsc->scalediag,&lsc->scalediag,PETSC_NULL);CHKERRQ(ierr); 113 } 114 ierr = PetscOptionsTail();CHKERRQ(ierr); 115 PetscFunctionReturn(0); 116 } 117 118 /*MC 119 PCLSC - Preconditioning for Schur complements, based on Least Squares Commutators 120 121 Options Database Key: 122 . -pc_lsc_scale_diag - Use the diagonal of A for scaling 123 124 Level: intermediate 125 126 Notes: 127 This preconditioner will normally be used with PCFieldSplit to precondition the Schur complement, but 128 it can be used for any Schur complement system. Consider the Schur complement 129 130 .vb 131 S = D - C inv(A) B 132 .ve 133 134 PCLSC currently doesn't do anything with D, so let's assume it is 0. The idea is that a good approximation to 135 inv(S) is given by 136 137 .vb 138 inv(CB) C A B inv(CB) 139 .ve 140 141 At some point, we'll be able to form the product CB for you, but for now the application has to provide it (this is 142 usually more efficient anyway). In the case of incompressible flow, CB is a Laplacian, call it L. The current 143 interface is to hang L and a preconditioning matrix Lp on the preconditioning matrix. 144 145 If you had called KSPSetOperators(ksp,S,Sp,flg), S should have type MATSCHURCOMPLEMENT and Sp can be any type you 146 like (PCLSC doesn't use it directly) but should have matrices composed with it, under the names "LSC_L" and "LSC_Lp". 147 For example, you might have setup code like this 148 149 .vb 150 PetscObjectCompose((PetscObject)Sp,"LSC_L",(PetscObject)L); 151 PetscObjectCompose((PetscObject)Sp,"LSC_Lp",(PetscObject)Lp); 152 .ve 153 154 And then your Jacobian assembly would look like 155 156 .vb 157 PetscObjectQuery((PetscObject)Sp,"LSC_L",(PetscObject*)&L); 158 PetscObjectQuery((PetscObject)Sp,"LSC_Lp",(PetscObject*)&Lp); 159 if (L) { assembly L } 160 if (Lp) { assemble Lp } 161 .ve 162 163 With this, you should be able to choose LSC preconditioning, using e.g. ML's algebraic multigrid to solve with L 164 165 .vb 166 -fieldsplit_1_pc_type lsc -fieldsplit_1_lsc_pc_type ml 167 .ve 168 169 Since we do not use the values in Sp, you can still put an assembled matrix there to use normal preconditioners. 170 171 Concepts: physics based preconditioners, block preconditioners 172 173 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, Block_Preconditioners, PCFIELDSPLIT, 174 PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(), PCFieldSplitSetType(), PCFieldSplitSetIS(), PCFieldSplitSchurPrecondition(), 175 MatCreateSchurComplement() 176 M*/ 177 178 EXTERN_C_BEGIN 179 #undef __FUNCT__ 180 #define __FUNCT__ "PCCreate_LSC" 181 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_LSC(PC pc) 182 { 183 PC_LSC *lsc; 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 ierr = PetscNewLog(pc,PC_LSC,&lsc);CHKERRQ(ierr); 188 pc->data = (void*)lsc; 189 190 pc->ops->apply = PCApply_LSC; 191 pc->ops->applytranspose = 0; 192 pc->ops->setup = PCSetUp_LSC; 193 pc->ops->destroy = PCDestroy_LSC; 194 pc->ops->setfromoptions = PCSetFromOptions_LSC; 195 pc->ops->view = 0; 196 pc->ops->applyrichardson = 0; 197 PetscFunctionReturn(0); 198 } 199 EXTERN_C_END 200