Actual source code: ex25.c
petsc-3.3-p7 2013-05-11
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: #include <petscdmda.h>
19: #include <petscksp.h>
21: static PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
22: static PetscErrorCode ComputeRHS(KSP,Vec,void*);
24: typedef struct {
25: PetscInt k;
26: PetscScalar e;
27: } AppCtx;
31: int main(int argc,char **argv)
32: {
34: KSP ksp;
35: DM da;
36: AppCtx user;
37: Mat A;
38: Vec b,b2;
39: Vec x;
40: PetscReal nrm;
42: PetscInitialize(&argc,&argv,(char *)0,help);
44: user.k = 1;
45: user.e = .99;
46: PetscOptionsGetInt(0,"-k",&user.k,0);
47: PetscOptionsGetScalar(0,"-e",&user.e,0);
49: KSPCreate(PETSC_COMM_WORLD,&ksp);
50: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-3,1,1,0,&da);
51: KSPSetDM(ksp,da);
52: KSPSetComputeRHS(ksp,ComputeRHS,&user);
53: KSPSetComputeOperators(ksp,ComputeMatrix,&user);
54: KSPSetFromOptions(ksp);
55: KSPSolve(ksp,PETSC_NULL,PETSC_NULL);
57: KSPGetOperators(ksp,&A,PETSC_NULL,PETSC_NULL);
58: KSPGetSolution(ksp,&x);
59: KSPGetRhs(ksp,&b);
60: VecDuplicate(b,&b2);
61: MatMult(A,x,b2);
62: VecAXPY(b2,-1.0,b);
63: VecNorm(b2,NORM_MAX,&nrm);
64: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",nrm);
66: VecDestroy(&b2);
67: KSPDestroy(&ksp);
68: DMDestroy(&da);
69: PetscFinalize();
71: return 0;
72: }
76: static PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
77: {
79: PetscInt mx,idx[2];
80: PetscScalar h,v[2];
81: DM da;
84: KSPGetDM(ksp,&da);
85: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
86: h = 1.0/((mx-1));
87: VecSet(b,h);
88: idx[0] = 0; idx[1] = mx -1;
89: v[0] = v[1] = 0.0;
90: VecSetValues(b,2,idx,v,INSERT_VALUES);
91: VecAssemblyBegin(b);
92: VecAssemblyEnd(b);
93: return(0);
94: }
98: static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,MatStructure *str,void *ctx)
99: {
100: AppCtx *user = (AppCtx*)ctx;
102: PetscInt i,mx,xm,xs;
103: PetscScalar v[3],h,xlow,xhigh;
104: MatStencil row,col[3];
105: DM da;
108: KSPGetDM(ksp,&da);
109: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
110: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
111: h = 1.0/(mx-1);
113: for(i=xs; i<xs+xm; i++){
114: row.i = i;
115: if (i==0 || i==mx-1){
116: v[0] = 2.0;
117: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
118: } else {
119: xlow = h*(PetscReal)i - .5*h;
120: xhigh = xlow + h;
121: v[0] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xlow))/h;col[0].i = i-1;
122: 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;
123: v[2] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[2].i = i+1;
124: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
125: }
126: }
127: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
128: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
129: return(0);
130: }