1 /* 2 Partial differential equation 3 4 d d u = 1, 0 < x < 1, 5 -- -- 6 dx dx 7 with boundary conditions 8 9 u = 0 for x = 0, x = 1 10 11 This uses multigrid to solve the linear system 12 13 Demonstrates how to build a DMSHELL for managing multigrid. The DMSHELL simply creates a 14 DMDA1d to construct all the needed PETSc objects. 15 16 */ 17 18 static char help[] = "Solves 1D constant coefficient Laplacian using DMSHELL and multigrid.\n\n"; 19 20 #include <petscdm.h> 21 #include <petscdmda.h> 22 #include <petscdmshell.h> 23 #include <petscksp.h> 24 25 static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *); 26 static PetscErrorCode ComputeRHS(KSP, Vec, void *); 27 static PetscErrorCode CreateMatrix(DM, Mat *); 28 static PetscErrorCode CreateGlobalVector(DM, Vec *); 29 static PetscErrorCode CreateLocalVector(DM, Vec *); 30 static PetscErrorCode Refine(DM, MPI_Comm, DM *); 31 static PetscErrorCode Coarsen(DM, MPI_Comm, DM *); 32 static PetscErrorCode CreateInterpolation(DM, DM, Mat *, Vec *); 33 static PetscErrorCode CreateRestriction(DM, DM, Mat *); 34 static PetscErrorCode Destroy(void *); 35 36 static PetscErrorCode MyDMShellCreate(MPI_Comm comm, DM da, DM *shell) 37 { 38 PetscFunctionBeginUser; 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 PetscFunctionReturn(PETSC_SUCCESS); 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 PetscFunctionBeginUser; 84 PetscCall(DMDestroy((DM *)&ctx)); 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 static PetscErrorCode CreateMatrix(DM shell, Mat *A) 89 { 90 DM da; 91 92 PetscFunctionBeginUser; 93 PetscCall(DMShellGetContext(shell, &da)); 94 PetscCall(DMCreateMatrix(da, A)); 95 PetscFunctionReturn(PETSC_SUCCESS); 96 } 97 98 static PetscErrorCode CreateInterpolation(DM dm1, DM dm2, Mat *mat, Vec *vec) 99 { 100 DM da1, da2; 101 102 PetscFunctionBeginUser; 103 PetscCall(DMShellGetContext(dm1, &da1)); 104 PetscCall(DMShellGetContext(dm2, &da2)); 105 PetscCall(DMCreateInterpolation(da1, da2, mat, vec)); 106 PetscFunctionReturn(PETSC_SUCCESS); 107 } 108 109 static PetscErrorCode CreateRestriction(DM dm1, DM dm2, Mat *mat) 110 { 111 DM da1, da2; 112 Mat tmat; 113 114 PetscFunctionBeginUser; 115 PetscCall(DMShellGetContext(dm1, &da1)); 116 PetscCall(DMShellGetContext(dm2, &da2)); 117 PetscCall(DMCreateInterpolation(da1, da2, &tmat, NULL)); 118 PetscCall(MatTranspose(tmat, MAT_INITIAL_MATRIX, mat)); 119 PetscCall(MatDestroy(&tmat)); 120 PetscFunctionReturn(PETSC_SUCCESS); 121 } 122 123 static PetscErrorCode CreateGlobalVector(DM shell, Vec *x) 124 { 125 DM da; 126 127 PetscFunctionBeginUser; 128 PetscCall(DMShellGetContext(shell, &da)); 129 PetscCall(DMCreateGlobalVector(da, x)); 130 PetscCall(VecSetDM(*x, shell)); 131 PetscFunctionReturn(PETSC_SUCCESS); 132 } 133 134 static PetscErrorCode CreateLocalVector(DM shell, Vec *x) 135 { 136 DM da; 137 138 PetscFunctionBeginUser; 139 PetscCall(DMShellGetContext(shell, &da)); 140 PetscCall(DMCreateLocalVector(da, x)); 141 PetscCall(VecSetDM(*x, shell)); 142 PetscFunctionReturn(PETSC_SUCCESS); 143 } 144 145 static PetscErrorCode Refine(DM shell, MPI_Comm comm, DM *dmnew) 146 { 147 DM da, dafine; 148 149 PetscFunctionBeginUser; 150 PetscCall(DMShellGetContext(shell, &da)); 151 PetscCall(DMRefine(da, comm, &dafine)); 152 PetscCall(MyDMShellCreate(PetscObjectComm((PetscObject)shell), dafine, dmnew)); 153 PetscFunctionReturn(PETSC_SUCCESS); 154 } 155 156 static PetscErrorCode Coarsen(DM shell, MPI_Comm comm, DM *dmnew) 157 { 158 DM da, dacoarse; 159 160 PetscFunctionBeginUser; 161 PetscCall(DMShellGetContext(shell, &da)); 162 PetscCall(DMCoarsen(da, comm, &dacoarse)); 163 PetscCall(MyDMShellCreate(PetscObjectComm((PetscObject)shell), dacoarse, dmnew)); 164 PetscFunctionReturn(PETSC_SUCCESS); 165 } 166 167 static PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx) 168 { 169 PetscInt mx, idx[2]; 170 PetscScalar h, v[2]; 171 DM da, shell; 172 173 PetscFunctionBeginUser; 174 PetscCall(KSPGetDM(ksp, &shell)); 175 PetscCall(DMShellGetContext(shell, &da)); 176 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 177 h = 1.0 / ((mx - 1)); 178 PetscCall(VecSet(b, h)); 179 idx[0] = 0; 180 idx[1] = mx - 1; 181 v[0] = v[1] = 0.0; 182 PetscCall(VecSetValues(b, 2, idx, v, INSERT_VALUES)); 183 PetscCall(VecAssemblyBegin(b)); 184 PetscCall(VecAssemblyEnd(b)); 185 PetscFunctionReturn(PETSC_SUCCESS); 186 } 187 188 static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx) 189 { 190 PetscInt i, mx, xm, xs; 191 PetscScalar v[3], h; 192 MatStencil row, col[3]; 193 DM da, shell; 194 195 PetscFunctionBeginUser; 196 PetscCall(KSPGetDM(ksp, &shell)); 197 PetscCall(DMShellGetContext(shell, &da)); 198 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 199 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 200 h = 1.0 / (mx - 1); 201 202 for (i = xs; i < xs + xm; i++) { 203 row.i = i; 204 if (i == 0 || i == mx - 1) { 205 v[0] = 2.0 / h; 206 PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES)); 207 } else { 208 v[0] = (-1.0) / h; 209 col[0].i = i - 1; 210 v[1] = (2.0) / h; 211 col[1].i = row.i; 212 v[2] = (-1.0) / h; 213 col[2].i = i + 1; 214 PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES)); 215 } 216 } 217 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 218 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 219 PetscFunctionReturn(PETSC_SUCCESS); 220 } 221 222 /*TEST 223 224 test: 225 nsize: 4 226 args: -ksp_monitor -pc_type mg -da_refine 3 227 228 TEST*/ 229