1 /* 2 Partial differential equation 3 4 d (1 + e*sine(2*pi*k*x)) 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 */ 14 15 static char help[] = "Solves 1D variable coefficient Laplacian using multigrid.\n\n"; 16 17 #include <petscdm.h> 18 #include <petscdmda.h> 19 #include <petscksp.h> 20 21 static PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *); 22 static PetscErrorCode ComputeRHS(KSP, Vec, void *); 23 24 typedef struct { 25 PetscInt k; 26 PetscScalar e; 27 } AppCtx; 28 29 int main(int argc, char **argv) 30 { 31 KSP ksp; 32 DM da; 33 AppCtx user; 34 Mat A; 35 Vec b, b2; 36 Vec x; 37 PetscReal nrm; 38 39 PetscFunctionBeginUser; 40 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 41 user.k = 1; 42 user.e = .99; 43 PetscCall(PetscOptionsGetInt(NULL, 0, "-k", &user.k, 0)); 44 PetscCall(PetscOptionsGetScalar(NULL, 0, "-e", &user.e, 0)); 45 46 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); 47 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 128, 1, 1, 0, &da)); 48 PetscCall(DMSetFromOptions(da)); 49 PetscCall(DMSetUp(da)); 50 PetscCall(KSPSetDM(ksp, da)); 51 PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, &user)); 52 PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, &user)); 53 PetscCall(KSPSetFromOptions(ksp)); 54 PetscCall(KSPSolve(ksp, NULL, NULL)); 55 56 PetscCall(KSPGetOperators(ksp, &A, NULL)); 57 PetscCall(KSPGetSolution(ksp, &x)); 58 PetscCall(KSPGetRhs(ksp, &b)); 59 PetscCall(VecDuplicate(b, &b2)); 60 PetscCall(MatMult(A, x, b2)); 61 PetscCall(VecAXPY(b2, -1.0, b)); 62 PetscCall(VecNorm(b2, NORM_MAX, &nrm)); 63 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)nrm)); 64 65 PetscCall(VecDestroy(&b2)); 66 PetscCall(KSPDestroy(&ksp)); 67 PetscCall(DMDestroy(&da)); 68 PetscCall(PetscFinalize()); 69 return 0; 70 } 71 72 static PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx) 73 { 74 PetscInt mx, idx[2]; 75 PetscScalar h, v[2]; 76 DM da; 77 78 PetscFunctionBeginUser; 79 PetscCall(KSPGetDM(ksp, &da)); 80 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 81 h = 1.0 / (mx - 1); 82 PetscCall(VecSet(b, h)); 83 idx[0] = 0; 84 idx[1] = mx - 1; 85 v[0] = v[1] = 0.0; 86 PetscCall(VecSetValues(b, 2, idx, v, INSERT_VALUES)); 87 PetscCall(VecAssemblyBegin(b)); 88 PetscCall(VecAssemblyEnd(b)); 89 PetscFunctionReturn(PETSC_SUCCESS); 90 } 91 92 static PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx) 93 { 94 AppCtx *user = (AppCtx *)ctx; 95 PetscInt i, mx, xm, xs; 96 PetscScalar v[3], h, xlow, xhigh; 97 MatStencil row, col[3]; 98 DM da; 99 100 PetscFunctionBeginUser; 101 PetscCall(KSPGetDM(ksp, &da)); 102 PetscCall(DMDAGetInfo(da, 0, &mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 103 PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0)); 104 h = 1.0 / (mx - 1); 105 106 for (i = xs; i < xs + xm; i++) { 107 row.i = i; 108 if (i == 0 || i == mx - 1) { 109 v[0] = 2.0 / h; 110 PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES)); 111 } else { 112 xlow = h * (PetscReal)i - .5 * h; 113 xhigh = xlow + h; 114 v[0] = (-1.0 - user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xlow)) / h; 115 col[0].i = i - 1; 116 v[1] = (2.0 + user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xlow) + user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xhigh)) / h; 117 col[1].i = row.i; 118 v[2] = (-1.0 - user->e * PetscSinScalar(2.0 * PETSC_PI * user->k * xhigh)) / h; 119 col[2].i = i + 1; 120 PetscCall(MatSetValuesStencil(jac, 1, &row, 3, col, v, INSERT_VALUES)); 121 } 122 } 123 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 124 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 125 PetscFunctionReturn(PETSC_SUCCESS); 126 } 127 128 /*TEST 129 130 test: 131 args: -pc_type mg -ksp_type fgmres -da_refine 2 -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -ksp_view -pc_mg_type full 132 requires: !single 133 134 test: 135 suffix: 2 136 nsize: 2 137 args: -pc_type mg -ksp_type fgmres -da_refine 2 -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -ksp_view -pc_mg_type full 138 requires: !single 139 140 TEST*/ 141