Actual source code: ex4.c
2: static char help[] = "Bilinear elements on the unit square for the Laplacian. Input arguments are:\n\
3: -m <size> : problem size\n\n";
5: #include <petscksp.h>
7: int FormElementStiffness(PetscReal H,PetscScalar *Ke)
8: {
9: Ke[0] = H/6.0; Ke[1] = -.125*H; Ke[2] = H/12.0; Ke[3] = -.125*H;
10: Ke[4] = -.125*H; Ke[5] = H/6.0; Ke[6] = -.125*H; Ke[7] = H/12.0;
11: Ke[8] = H/12.0; Ke[9] = -.125*H; Ke[10] = H/6.0; Ke[11] = -.125*H;
12: Ke[12] = -.125*H; Ke[13] = H/12.0; Ke[14] = -.125*H; Ke[15] = H/6.0;
13: return 0;
14: }
16: int FormElementRhs(PetscReal x,PetscReal y,PetscReal H,PetscScalar *r)
17: {
18: r[0] = 0.; r[1] = 0.; r[2] = 0.; r[3] = 0.0;
19: return 0;
20: }
22: /* Note: this code is for testing purposes only. The assembly process is not scalable */
23: int main(int argc,char **args)
24: {
25: Mat C;
26: PetscInt i,m = 2,N,M,its,idx[4],count,*rows;
27: PetscScalar val,Ke[16],r[4];
28: PetscReal x,y,h,norm;
29: Vec u,ustar,b;
30: KSP ksp;
31: PetscMPIInt rank;
32: PetscBool usezerorows = PETSC_TRUE;
34: PetscInitialize(&argc,&args,(char*)0,help);
35: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
36: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
37: PetscOptionsGetBool(NULL,NULL,"-usezerorows",&usezerorows,NULL);
38: N = (m+1)*(m+1); /* dimension of matrix */
39: M = m*m; /* number of elements */
40: h = 1.0/m; /* mesh width */
42: /* create stiffness matrix */
43: MatCreate(PETSC_COMM_WORLD,&C);
44: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);
45: MatSetFromOptions(C);
46: MatMPIAIJSetPreallocation(C,9,NULL,9,NULL);
47: MatSeqAIJSetPreallocation(C,9,NULL);
48: #if defined(PETSC_HAVE_HYPRE)
49: MatHYPRESetPreallocation(C,9,NULL,9,NULL);
50: #endif
52: /* forms the element stiffness for the Laplacian */
53: FormElementStiffness(h*h,Ke);
55: /* assemble the matrix: only process 0 adds the values, not scalable */
56: if (!rank) {
57: for (i=0; i<M; i++) {
58: /* node numbers for the four corners of element */
59: idx[0] = (m+1)*(i/m) + (i % m);
60: idx[1] = idx[0] + 1;
61: idx[2] = idx[1] + m + 1;
62: idx[3] = idx[2] - 1;
63: if (i == M-1 && !usezerorows) { /* If MatZeroRows not supported -> make it non-singular */
64: for (PetscInt ii = 0; ii < 4; ii++) {
65: Ke[ 2*4 + ii] = 0.0;
66: Ke[ii*4 + 2] = 0.0;
67: }
68: Ke[10] = 1.0;
69: }
70: MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
71: }
72: }
73: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
76: /* create right hand side and solution */
77: MatCreateVecs(C,&u,&b);
78: VecDuplicate(u,&ustar);
79: VecSet(u,0.0);
80: VecSet(b,0.0);
82: /* assemble the right hand side: only process 0 adds the values, not scalable */
83: if (!rank) {
84: for (i=0; i<M; i++) {
85: /* location of lower left corner of element */
86: x = h*(i%m);
87: y = h*(i/m);
88: /* node numbers for the four corners of element */
89: idx[0] = (m+1)*(i/m) + (i%m);
90: idx[1] = idx[0]+1;
91: idx[2] = idx[1]+m+1;
92: idx[3] = idx[2]-1;
93: FormElementRhs(x,y,h*h,r);
94: VecSetValues(b,4,idx,r,ADD_VALUES);
95: }
96: }
97: VecAssemblyBegin(b);
98: VecAssemblyEnd(b);
100: /* modify matrix and rhs for Dirichlet boundary conditions */
101: if (!rank) {
102: PetscMalloc1(4*m+1,&rows);
103: for (i=0; i<m+1; i++) {
104: rows[i] = i; /* bottom */
105: rows[3*m-1+i] = m*(m+1)+i; /* top */
106: }
107: count = m+1; /* left side */
108: for (i=m+1; i<m*(m+1); i+=m+1) rows[count++] = i;
110: count = 2*m; /* right side */
111: for (i=2*m+1; i<m*(m+1); i+=m+1) rows[count++] = i;
113: for (i=0; i<4*m; i++) {
114: val = h*(rows[i]/(m+1));
115: VecSetValues(u,1,&rows[i],&val,INSERT_VALUES);
116: VecSetValues(b,1,&rows[i],&val,INSERT_VALUES);
117: }
118: if (usezerorows) MatZeroRows(C,4*m,rows,1.0,0,0);
119: PetscFree(rows);
120: } else {
121: if (usezerorows) MatZeroRows(C,0,NULL,1.0,0,0);
122: }
123: VecAssemblyBegin(u);
124: VecAssemblyEnd(u);
125: VecAssemblyBegin(b);
126: VecAssemblyEnd(b);
128: if (!usezerorows) {
129: VecSet(ustar,1.0);
130: MatMult(C,ustar,b);
131: }
133: /* solve linear system */
134: KSPCreate(PETSC_COMM_WORLD,&ksp);
135: KSPSetOperators(ksp,C,C);
136: KSPSetFromOptions(ksp);
137: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
138: KSPSolve(ksp,b,u);
140: /* check error */
141: if (usezerorows) {
142: if (!rank) {
143: for (i=0; i<N; i++) {
144: val = h*(i/(m+1));
145: VecSetValues(ustar,1,&i,&val,INSERT_VALUES);
146: }
147: }
148: VecAssemblyBegin(ustar);
149: VecAssemblyEnd(ustar);
150: }
152: VecAXPY(u,-1.0,ustar);
153: VecNorm(u,NORM_2,&norm);
154: KSPGetIterationNumber(ksp,&its);
155: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)(norm*h),its);
157: KSPDestroy(&ksp);
158: VecDestroy(&ustar);
159: VecDestroy(&u);
160: VecDestroy(&b);
161: MatDestroy(&C);
162: PetscFinalize();
163: return 0;
164: }
166: /*TEST
168: test:
169: args: -ksp_monitor_short -m 5 -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
171: test:
172: suffix: 3
173: args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
175: test:
176: suffix: 5
177: args: -pc_type eisenstat -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
179: test:
180: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
181: suffix: hypre_device_none
182: output_file: output/ex4_hypre_none.out
183: nsize: {{1 2}}
184: args: -usezerorows 0 -mat_type hypre -pc_type none -m 5
186: test:
187: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
188: suffix: hypre_device_amg
189: nsize: {{1 2}}
190: args: -usezerorows 0 -mat_type hypre -pc_type hypre -m 25 -ksp_type cg -ksp_norm_type natural
192: TEST*/