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