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 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; idx[1] = mx -1; 84 v[0] = v[1] = 0.0; 85 PetscCall(VecSetValues(b,2,idx,v,INSERT_VALUES)); 86 PetscCall(VecAssemblyBegin(b)); 87 PetscCall(VecAssemblyEnd(b)); 88 PetscFunctionReturn(0); 89 } 90 91 static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx) 92 { 93 AppCtx *user = (AppCtx*)ctx; 94 PetscInt i,mx,xm,xs; 95 PetscScalar v[3],h,xlow,xhigh; 96 MatStencil row,col[3]; 97 DM da; 98 99 PetscFunctionBeginUser; 100 PetscCall(KSPGetDM(ksp,&da)); 101 PetscCall(DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0)); 102 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 103 h = 1.0/(mx-1); 104 105 for (i=xs; i<xs+xm; i++) { 106 row.i = i; 107 if (i==0 || i==mx-1) { 108 v[0] = 2.0/h; 109 PetscCall(MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES)); 110 } else { 111 xlow = h*(PetscReal)i - .5*h; 112 xhigh = xlow + h; 113 v[0] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xlow))/h;col[0].i = i-1; 114 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; 115 v[2] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[2].i = i+1; 116 PetscCall(MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES)); 117 } 118 } 119 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 120 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 121 PetscFunctionReturn(0); 122 } 123 124 /*TEST 125 126 test: 127 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 128 requires: !single 129 130 test: 131 suffix: 2 132 nsize: 2 133 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 134 requires: !single 135 136 TEST*/ 137