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