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