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