1 2 /* 3 Partial differential equation 4 5 d d u = 1, 0 < x < 1, 6 -- -- 7 dx dx 8 with boundary conditions 9 10 u = 0 for x = 0, x = 1 11 12 This uses multigrid to solve the linear system 13 14 Demonstrates how to build a DMSHELL for managing multigrid. The DMSHELL simply creates a 15 DMDA1d to construct all the needed PETSc objects. 16 17 */ 18 19 static char help[] = "Solves 1D constant coefficient Laplacian using DMSHELL and multigrid.\n\n"; 20 21 #include <petscdm.h> 22 #include <petscdmda.h> 23 #include <petscdmshell.h> 24 #include <petscksp.h> 25 26 static PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*); 27 static PetscErrorCode ComputeRHS(KSP,Vec,void*); 28 static PetscErrorCode CreateMatrix(DM,Mat*); 29 static PetscErrorCode CreateGlobalVector(DM,Vec*); 30 static PetscErrorCode CreateLocalVector(DM,Vec*); 31 static PetscErrorCode Refine(DM,MPI_Comm,DM*); 32 static PetscErrorCode Coarsen(DM,MPI_Comm,DM*); 33 static PetscErrorCode CreateInterpolation(DM,DM,Mat*,Vec*); 34 static PetscErrorCode CreateRestriction(DM,DM,Mat*); 35 36 static PetscErrorCode MyDMShellCreate(MPI_Comm comm,DM da,DM *shell) 37 { 38 PetscErrorCode ierr; 39 40 ierr = DMShellCreate(comm,shell);CHKERRQ(ierr); 41 ierr = DMShellSetContext(*shell,da);CHKERRQ(ierr); 42 ierr = DMShellSetCreateMatrix(*shell,CreateMatrix);CHKERRQ(ierr); 43 ierr = DMShellSetCreateGlobalVector(*shell,CreateGlobalVector);CHKERRQ(ierr); 44 ierr = DMShellSetCreateLocalVector(*shell,CreateLocalVector);CHKERRQ(ierr); 45 ierr = DMShellSetRefine(*shell,Refine);CHKERRQ(ierr); 46 ierr = DMShellSetCoarsen(*shell,Coarsen);CHKERRQ(ierr); 47 ierr = DMShellSetCreateInterpolation(*shell,CreateInterpolation);CHKERRQ(ierr); 48 ierr = DMShellSetCreateRestriction(*shell,CreateRestriction);CHKERRQ(ierr); 49 return 0; 50 } 51 52 int main(int argc,char **argv) 53 { 54 PetscErrorCode ierr; 55 KSP ksp; 56 DM da,shell; 57 PetscInt levels; 58 59 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 60 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); 61 ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,129,1,1,0,&da);CHKERRQ(ierr); 62 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 63 ierr = DMSetUp(da);CHKERRQ(ierr); 64 ierr = MyDMShellCreate(PETSC_COMM_WORLD,da,&shell);CHKERRQ(ierr); 65 /* these two lines are not needed but allow PCMG to automatically know how many multigrid levels the user wants */ 66 ierr = DMGetRefineLevel(da,&levels);CHKERRQ(ierr); 67 ierr = DMSetRefineLevel(shell,levels);CHKERRQ(ierr); 68 69 ierr = KSPSetDM(ksp,shell);CHKERRQ(ierr); 70 ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr); 71 ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);CHKERRQ(ierr); 72 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); 73 ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr); 74 75 ierr = KSPDestroy(&ksp);CHKERRQ(ierr); 76 ierr = DMDestroy(&shell);CHKERRQ(ierr); 77 ierr = DMDestroy(&da);CHKERRQ(ierr); 78 ierr = PetscFinalize(); 79 return ierr; 80 } 81 82 static PetscErrorCode CreateMatrix(DM shell,Mat *A) 83 { 84 PetscErrorCode ierr; 85 DM da; 86 87 ierr = DMShellGetContext(shell,(void**)&da);CHKERRQ(ierr); 88 ierr = DMCreateMatrix(da,A);CHKERRQ(ierr); 89 return 0; 90 } 91 92 static PetscErrorCode CreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec) 93 { 94 DM da1,da2; 95 PetscErrorCode ierr; 96 97 ierr = DMShellGetContext(dm1,(void**)&da1);CHKERRQ(ierr); 98 ierr = DMShellGetContext(dm2,(void**)&da2);CHKERRQ(ierr); 99 ierr = DMCreateInterpolation(da1,da2,mat,vec);CHKERRQ(ierr); 100 return 0; 101 } 102 103 static PetscErrorCode CreateRestriction(DM dm1,DM dm2,Mat *mat) 104 { 105 DM da1,da2; 106 PetscErrorCode ierr; 107 Mat tmat; 108 109 ierr = DMShellGetContext(dm1,(void**)&da1);CHKERRQ(ierr); 110 ierr = DMShellGetContext(dm2,(void**)&da2);CHKERRQ(ierr); 111 ierr = DMCreateInterpolation(da1,da2,&tmat,NULL);CHKERRQ(ierr); 112 ierr = MatTranspose(tmat,MAT_INITIAL_MATRIX,mat);CHKERRQ(ierr); 113 ierr = MatDestroy(&tmat);CHKERRQ(ierr); 114 return 0; 115 } 116 117 static PetscErrorCode CreateGlobalVector(DM shell,Vec *x) 118 { 119 PetscErrorCode ierr; 120 DM da; 121 122 ierr = DMShellGetContext(shell,(void**)&da);CHKERRQ(ierr); 123 ierr = DMCreateGlobalVector(da,x);CHKERRQ(ierr); 124 ierr = VecSetDM(*x,shell);CHKERRQ(ierr); 125 return 0; 126 } 127 128 static PetscErrorCode CreateLocalVector(DM shell,Vec *x) 129 { 130 PetscErrorCode ierr; 131 DM da; 132 133 ierr = DMShellGetContext(shell,(void**)&da);CHKERRQ(ierr); 134 ierr = DMCreateLocalVector(da,x);CHKERRQ(ierr); 135 ierr = VecSetDM(*x,shell);CHKERRQ(ierr); 136 return 0; 137 } 138 139 static PetscErrorCode Refine(DM shell,MPI_Comm comm,DM *dmnew) 140 { 141 PetscErrorCode ierr; 142 DM da,dafine; 143 144 ierr = DMShellGetContext(shell,(void**)&da);CHKERRQ(ierr); 145 ierr = DMRefine(da,comm,&dafine);CHKERRQ(ierr); 146 ierr = MyDMShellCreate(PetscObjectComm((PetscObject)shell),dafine,dmnew);CHKERRQ(ierr); 147 return 0; 148 } 149 150 static PetscErrorCode Coarsen(DM shell,MPI_Comm comm,DM *dmnew) 151 { 152 PetscErrorCode ierr; 153 DM da,dacoarse; 154 155 ierr = DMShellGetContext(shell,(void**)&da);CHKERRQ(ierr); 156 ierr = DMCoarsen(da,comm,&dacoarse);CHKERRQ(ierr); 157 ierr = MyDMShellCreate(PetscObjectComm((PetscObject)shell),dacoarse,dmnew);CHKERRQ(ierr); 158 /* discard an "extra" reference count to dacoarse */ 159 ierr = DMDestroy(&dacoarse);CHKERRQ(ierr); 160 return 0; 161 } 162 163 static PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx) 164 { 165 PetscErrorCode ierr; 166 PetscInt mx,idx[2]; 167 PetscScalar h,v[2]; 168 DM da,shell; 169 170 PetscFunctionBeginUser; 171 ierr = KSPGetDM(ksp,&shell);CHKERRQ(ierr); 172 ierr = DMShellGetContext(shell,(void**)&da);CHKERRQ(ierr); 173 ierr = DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 174 h = 1.0/((mx-1)); 175 ierr = VecSet(b,h);CHKERRQ(ierr); 176 idx[0] = 0; idx[1] = mx -1; 177 v[0] = v[1] = 0.0; 178 ierr = VecSetValues(b,2,idx,v,INSERT_VALUES);CHKERRQ(ierr); 179 ierr = VecAssemblyBegin(b);CHKERRQ(ierr); 180 ierr = VecAssemblyEnd(b);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx) 185 { 186 PetscErrorCode ierr; 187 PetscInt i,mx,xm,xs; 188 PetscScalar v[3],h; 189 MatStencil row,col[3]; 190 DM da,shell; 191 192 PetscFunctionBeginUser; 193 ierr = KSPGetDM(ksp,&shell);CHKERRQ(ierr); 194 ierr = DMShellGetContext(shell,(void**)&da);CHKERRQ(ierr); 195 ierr = DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 196 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 197 h = 1.0/(mx-1); 198 199 for (i=xs; i<xs+xm; i++) { 200 row.i = i; 201 if (i==0 || i==mx-1) { 202 v[0] = 2.0/h; 203 ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr); 204 } else { 205 v[0] = (-1.0)/h;col[0].i = i-1; 206 v[1] = (2.0)/h;col[1].i = row.i; 207 v[2] = (-1.0)/h;col[2].i = i+1; 208 ierr = MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);CHKERRQ(ierr); 209 } 210 } 211 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 212 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 217 /*TEST 218 219 test: 220 nsize: 4 221 args: -ksp_monitor -pc_type mg -da_refine 3 222 223 TEST*/ 224