Actual source code: ex25.c
2: /*
3: Partial differential equation
5: d (1 + e*sine(2*pi*k*x)) d u = 1, 0 < x < 1,
6: -- ---
7: dx dx
8: with boundary conditions
10: u = 0 for x = 0, x = 1
12: This uses multigrid to solve the linear system
14: */
16: static char help[] = "Solves 1D variable coefficient Laplacian using multigrid.\n\n";
18: /*T
20: T*/
22: #include <petscdm.h>
23: #include <petscdmda.h>
24: #include <petscksp.h>
26: static PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
27: static PetscErrorCode ComputeRHS(KSP,Vec,void*);
29: typedef struct {
30: PetscInt k;
31: PetscScalar e;
32: } AppCtx;
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;
44: PetscInitialize(&argc,&argv,(char*)0,help);
45: user.k = 1;
46: user.e = .99;
47: PetscOptionsGetInt(NULL,0,"-k",&user.k,0);
48: PetscOptionsGetScalar(NULL,0,"-e",&user.e,0);
50: KSPCreate(PETSC_COMM_WORLD,&ksp);
51: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,128,1,1,0,&da);
52: DMSetFromOptions(da);
53: DMSetUp(da);
54: KSPSetDM(ksp,da);
55: KSPSetComputeRHS(ksp,ComputeRHS,&user);
56: KSPSetComputeOperators(ksp,ComputeMatrix,&user);
57: KSPSetFromOptions(ksp);
58: KSPSolve(ksp,NULL,NULL);
60: KSPGetOperators(ksp,&A,NULL);
61: KSPGetSolution(ksp,&x);
62: KSPGetRhs(ksp,&b);
63: VecDuplicate(b,&b2);
64: MatMult(A,x,b2);
65: VecAXPY(b2,-1.0,b);
66: VecNorm(b2,NORM_MAX,&nrm);
67: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)nrm);
69: VecDestroy(&b2);
70: KSPDestroy(&ksp);
71: DMDestroy(&da);
72: PetscFinalize();
73: return 0;
74: }
76: static PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
77: {
78: PetscInt mx,idx[2];
79: PetscScalar h,v[2];
80: DM da;
83: KSPGetDM(ksp,&da);
84: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
85: h = 1.0/((mx-1));
86: VecSet(b,h);
87: idx[0] = 0; idx[1] = mx -1;
88: v[0] = v[1] = 0.0;
89: VecSetValues(b,2,idx,v,INSERT_VALUES);
90: VecAssemblyBegin(b);
91: VecAssemblyEnd(b);
92: return 0;
93: }
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;
104: KSPGetDM(ksp,&da);
105: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
106: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
107: h = 1.0/(mx-1);
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: 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: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
121: }
122: }
123: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
124: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
125: return 0;
126: }
128: /*TEST
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
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
140: TEST*/