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 static PetscErrorCode Destroy(void *); 36 37 static PetscErrorCode MyDMShellCreate(MPI_Comm comm, DM da, DM *shell) 38 { 39 PetscCall(DMShellCreate(comm, shell)); 40 PetscCall(DMShellSetContext(*shell, da)); 41 PetscCall(DMShellSetCreateMatrix(*shell, CreateMatrix)); 42 PetscCall(DMShellSetCreateGlobalVector(*shell, CreateGlobalVector)); 43 PetscCall(DMShellSetCreateLocalVector(*shell, CreateLocalVector)); 44 PetscCall(DMShellSetRefine(*shell, Refine)); 45 PetscCall(DMShellSetCoarsen(*shell, Coarsen)); 46 PetscCall(DMShellSetCreateInterpolation(*shell, CreateInterpolation)); 47 PetscCall(DMShellSetCreateRestriction(*shell, CreateRestriction)); 48 PetscCall(DMShellSetDestroyContext(*shell, Destroy)); 49 return 0; 50 } 51 52 int main(int argc, char **argv) 53 { 54 KSP ksp; 55 DM da, shell; 56 PetscInt levels; 57 58 PetscFunctionBeginUser; 59 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 60 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 61 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 129, 1, 1, 0, &da)); 62 PetscCall(DMSetFromOptions(da)); 63 PetscCall(DMSetUp(da)); 64 PetscCall(MyDMShellCreate(PETSC_COMM_WORLD, da, &shell)); 65 /* these two lines are not needed but allow PCMG to automatically know how many multigrid levels the user wants */ 66 PetscCall(DMGetRefineLevel(da, &levels)); 67 PetscCall(DMSetRefineLevel(shell, levels)); 68 69 PetscCall(KSPSetDM(ksp, shell)); 70 PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, NULL)); 71 PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, NULL)); 72 PetscCall(KSPSetFromOptions(ksp)); 73 PetscCall(KSPSolve(ksp, NULL, NULL)); 74 75 PetscCall(KSPDestroy(&ksp)); 76 PetscCall(DMDestroy(&shell)); 77 PetscCall(PetscFinalize()); 78 return 0; 79 } 80 81 static PetscErrorCode Destroy(void *ctx) 82 { 83 PetscCall(DMDestroy((DM *)&ctx)); 84 return 0; 85 } 86 87 static PetscErrorCode CreateMatrix(DM shell, Mat *A) 88 { 89 DM da; 90 91 PetscCall(DMShellGetContext(shell, &da)); 92 PetscCall(DMCreateMatrix(da, A)); 93 return 0; 94 } 95 96 static PetscErrorCode CreateInterpolation(DM dm1, DM dm2, Mat *mat, Vec *vec) 97 { 98 DM da1, da2; 99 100 PetscCall(DMShellGetContext(dm1, &da1)); 101 PetscCall(DMShellGetContext(dm2, &da2)); 102 PetscCall(DMCreateInterpolation(da1, da2, mat, vec)); 103 return 0; 104 } 105 106 static PetscErrorCode CreateRestriction(DM dm1, DM dm2, Mat *mat) 107 { 108 DM da1, da2; 109 Mat tmat; 110 111 PetscCall(DMShellGetContext(dm1, &da1)); 112 PetscCall(DMShellGetContext(dm2, &da2)); 113 PetscCall(DMCreateInterpolation(da1, da2, &tmat, NULL)); 114 PetscCall(MatTranspose(tmat, MAT_INITIAL_MATRIX, mat)); 115 PetscCall(MatDestroy(&tmat)); 116 return 0; 117 } 118 119 static PetscErrorCode CreateGlobalVector(DM shell, Vec *x) 120 { 121 DM da; 122 123 PetscCall(DMShellGetContext(shell, &da)); 124 PetscCall(DMCreateGlobalVector(da, x)); 125 PetscCall(VecSetDM(*x, shell)); 126 return 0; 127 } 128 129 static PetscErrorCode CreateLocalVector(DM shell, Vec *x) 130 { 131 DM da; 132 133 PetscCall(DMShellGetContext(shell, &da)); 134 PetscCall(DMCreateLocalVector(da, x)); 135 PetscCall(VecSetDM(*x, shell)); 136 return 0; 137 } 138 139 static PetscErrorCode Refine(DM shell, MPI_Comm comm, DM *dmnew) 140 { 141 DM da, dafine; 142 143 PetscCall(DMShellGetContext(shell, &da)); 144 PetscCall(DMRefine(da, comm, &dafine)); 145 PetscCall(MyDMShellCreate(PetscObjectComm((PetscObject)shell), dafine, dmnew)); 146 return 0; 147 } 148 149 static PetscErrorCode Coarsen(DM shell, MPI_Comm comm, DM *dmnew) 150 { 151 DM da, dacoarse; 152 153 PetscCall(DMShellGetContext(shell, &da)); 154 PetscCall(DMCoarsen(da, comm, &dacoarse)); 155 PetscCall(MyDMShellCreate(PetscObjectComm((PetscObject)shell), dacoarse, dmnew)); 156 return 0; 157 } 158 159 static PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx) 160 { 161 PetscInt mx, idx[2]; 162 PetscScalar h, v[2]; 163 DM da, shell; 164 165 PetscFunctionBeginUser; 166 PetscCall(KSPGetDM(ksp, &shell)); 167 PetscCall(DMShellGetContext(shell, &da)); 168 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 169 h = 1.0 / ((mx - 1)); 170 PetscCall(VecSet(b, h)); 171 idx[0] = 0; 172 idx[1] = mx - 1; 173 v[0] = v[1] = 0.0; 174 PetscCall(VecSetValues(b, 2, idx, v, INSERT_VALUES)); 175 PetscCall(VecAssemblyBegin(b)); 176 PetscCall(VecAssemblyEnd(b)); 177 PetscFunctionReturn(0); 178 } 179 180 static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx) 181 { 182 PetscInt i, mx, xm, xs; 183 PetscScalar v[3], h; 184 MatStencil row, col[3]; 185 DM da, shell; 186 187 PetscFunctionBeginUser; 188 PetscCall(KSPGetDM(ksp, &shell)); 189 PetscCall(DMShellGetContext(shell, &da)); 190 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 191 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 192 h = 1.0 / (mx - 1); 193 194 for (i = xs; i < xs + xm; i++) { 195 row.i = i; 196 if (i == 0 || i == mx - 1) { 197 v[0] = 2.0 / h; 198 PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES)); 199 } else { 200 v[0] = (-1.0) / h; 201 col[0].i = i - 1; 202 v[1] = (2.0) / h; 203 col[1].i = row.i; 204 v[2] = (-1.0) / h; 205 col[2].i = i + 1; 206 PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES)); 207 } 208 } 209 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 210 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 211 PetscFunctionReturn(0); 212 } 213 214 /*TEST 215 216 test: 217 nsize: 4 218 args: -ksp_monitor -pc_type mg -da_refine 3 219 220 TEST*/ 221