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; 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(0); 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;col[0].i = i-1; 115 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; 116 v[2] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[2].i = i+1; 117 PetscCall(MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES)); 118 } 119 } 120 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 121 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 122 PetscFunctionReturn(0); 123 } 124 125 /*TEST 126 127 test: 128 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 129 requires: !single 130 131 test: 132 suffix: 2 133 nsize: 2 134 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 135 requires: !single 136 137 TEST*/ 138