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;
27: PetscInt i,m = 2,N,M,its,idx[4],count,*rows;
28: PetscScalar val,Ke[16],r[4];
29: PetscReal x,y,h,norm;
30: Vec u,ustar,b;
31: KSP ksp;
32: PetscMPIInt rank;
33: PetscBool usezerorows = PETSC_TRUE;
35: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
36: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
37: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
38: PetscOptionsGetBool(NULL,NULL,"-usezerorows",&usezerorows,NULL);
39: N = (m+1)*(m+1); /* dimension of matrix */
40: M = m*m; /* number of elements */
41: h = 1.0/m; /* mesh width */
43: /* create stiffness matrix */
44: MatCreate(PETSC_COMM_WORLD,&C);
45: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);
46: MatSetFromOptions(C);
47: MatMPIAIJSetPreallocation(C,9,NULL,9,NULL);
48: MatSeqAIJSetPreallocation(C,9,NULL);
49: #if defined(PETSC_HAVE_HYPRE)
50: MatHYPRESetPreallocation(C,9,NULL,9,NULL);
51: #endif
53: /* forms the element stiffness for the Laplacian */
54: FormElementStiffness(h*h,Ke);
56: /* assemble the matrix: only process 0 adds the values, not scalable */
57: if (!rank) {
58: for (i=0; i<M; i++) {
59: /* node numbers for the four corners of element */
60: idx[0] = (m+1)*(i/m) + (i % m);
61: idx[1] = idx[0] + 1;
62: idx[2] = idx[1] + m + 1;
63: idx[3] = idx[2] - 1;
64: if (i == M-1 && !usezerorows) { /* If MatZeroRows not supported -> make it non-singular */
65: for (PetscInt ii = 0; ii < 4; ii++) {
66: Ke[ 2*4 + ii] = 0.0;
67: Ke[ii*4 + 2] = 0.0;
68: }
69: Ke[10] = 1.0;
70: }
71: MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);
72: }
73: }
74: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
77: /* create right hand side and solution */
78: MatCreateVecs(C,&u,&b);
79: VecDuplicate(u,&ustar);
80: VecSet(u,0.0);
81: VecSet(b,0.0);
83: /* assemble the right hand side: only process 0 adds the values, not scalable */
84: if (!rank) {
85: for (i=0; i<M; i++) {
86: /* location of lower left corner of element */
87: x = h*(i%m);
88: y = h*(i/m);
89: /* node numbers for the four corners of element */
90: idx[0] = (m+1)*(i/m) + (i%m);
91: idx[1] = idx[0]+1;
92: idx[2] = idx[1]+m+1;
93: idx[3] = idx[2]-1;
94: FormElementRhs(x,y,h*h,r);
95: VecSetValues(b,4,idx,r,ADD_VALUES);
96: }
97: }
98: VecAssemblyBegin(b);
99: VecAssemblyEnd(b);
101: /* modify matrix and rhs for Dirichlet boundary conditions */
102: if (!rank) {
103: PetscMalloc1(4*m+1,&rows);
104: for (i=0; i<m+1; i++) {
105: rows[i] = i; /* bottom */
106: rows[3*m-1+i] = m*(m+1)+i; /* top */
107: }
108: count = m+1; /* left side */
109: for (i=m+1; i<m*(m+1); i+=m+1) rows[count++] = i;
111: count = 2*m; /* right side */
112: for (i=2*m+1; i<m*(m+1); i+=m+1) rows[count++] = i;
114: for (i=0; i<4*m; i++) {
115: val = h*(rows[i]/(m+1));
116: VecSetValues(u,1,&rows[i],&val,INSERT_VALUES);
117: VecSetValues(b,1,&rows[i],&val,INSERT_VALUES);
118: }
119: if (usezerorows) { MatZeroRows(C,4*m,rows,1.0,0,0); }
120: PetscFree(rows);
121: } else {
122: if (usezerorows) { MatZeroRows(C,0,NULL,1.0,0,0); }
123: }
124: VecAssemblyBegin(u);
125: VecAssemblyEnd(u);
126: VecAssemblyBegin(b);
127: VecAssemblyEnd(b);
129: if (!usezerorows) {
130: VecSet(ustar,1.0);
131: MatMult(C,ustar,b);
132: }
134: /* solve linear system */
135: KSPCreate(PETSC_COMM_WORLD,&ksp);
136: KSPSetOperators(ksp,C,C);
137: KSPSetFromOptions(ksp);
138: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
139: KSPSolve(ksp,b,u);
141: /* check error */
142: if (usezerorows) {
143: if (!rank) {
144: for (i=0; i<N; i++) {
145: val = h*(i/(m+1));
146: VecSetValues(ustar,1,&i,&val,INSERT_VALUES);
147: }
148: }
149: VecAssemblyBegin(ustar);
150: VecAssemblyEnd(ustar);
151: }
153: VecAXPY(u,-1.0,ustar);
154: VecNorm(u,NORM_2,&norm);
155: KSPGetIterationNumber(ksp,&its);
156: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)(norm*h),its);
158: KSPDestroy(&ksp);
159: VecDestroy(&ustar);
160: VecDestroy(&u);
161: VecDestroy(&b);
162: MatDestroy(&C);
163: PetscFinalize();
164: return ierr;
165: }
167: /*TEST
169: test:
170: args: -ksp_monitor_short -m 5 -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
172: test:
173: suffix: 3
174: args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
176: test:
177: suffix: 5
178: args: -pc_type eisenstat -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always
180: test:
181: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
182: suffix: hypre_device_none
183: output_file: output/ex4_hypre_none.out
184: nsize: {{1 2}}
185: args: -usezerorows 0 -mat_type hypre -pc_type none -m 5
187: test:
188: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
189: suffix: hypre_device_amg
190: nsize: {{1 2}}
191: args: -usezerorows 0 -mat_type hypre -pc_type hypre -m 25 -ksp_type cg -ksp_norm_type natural
193: TEST*/