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;
 10:   Ke[1]  = -.125 * H;
 11:   Ke[2]  = H / 12.0;
 12:   Ke[3]  = -.125 * H;
 13:   Ke[4]  = -.125 * H;
 14:   Ke[5]  = H / 6.0;
 15:   Ke[6]  = -.125 * H;
 16:   Ke[7]  = H / 12.0;
 17:   Ke[8]  = H / 12.0;
 18:   Ke[9]  = -.125 * H;
 19:   Ke[10] = H / 6.0;
 20:   Ke[11] = -.125 * H;
 21:   Ke[12] = -.125 * H;
 22:   Ke[13] = H / 12.0;
 23:   Ke[14] = -.125 * H;
 24:   Ke[15] = H / 6.0;
 25:   return 0;
 26: }

 28: int FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
 29: {
 30:   r[0] = 0.;
 31:   r[1] = 0.;
 32:   r[2] = 0.;
 33:   r[3] = 0.0;
 34:   return 0;
 35: }

 37: /* Note: this code is for testing purposes only. The assembly process is not scalable */
 38: int main(int argc, char **args)
 39: {
 40:   Mat         C;
 41:   PetscInt    i, m = 2, N, M, its, idx[4], count, *rows;
 42:   PetscScalar val, Ke[16], r[4];
 43:   PetscReal   x, y, h, norm;
 44:   Vec         u, ustar, b;
 45:   KSP         ksp;
 46:   PetscMPIInt rank;
 47:   PetscBool   usezerorows = PETSC_TRUE;

 50:   PetscInitialize(&argc, &args, (char *)0, help);
 51:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 52:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 53:   PetscOptionsGetBool(NULL, NULL, "-usezerorows", &usezerorows, NULL);
 54:   N = (m + 1) * (m + 1); /* dimension of matrix */
 55:   M = m * m;             /* number of elements */
 56:   h = 1.0 / m;           /* mesh width */

 58:   /* create stiffness matrix */
 59:   MatCreate(PETSC_COMM_WORLD, &C);
 60:   MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N);
 61:   MatSetFromOptions(C);
 62:   MatMPIAIJSetPreallocation(C, 9, NULL, 9, NULL);
 63:   MatSeqAIJSetPreallocation(C, 9, NULL);
 64: #if defined(PETSC_HAVE_HYPRE)
 65:   MatHYPRESetPreallocation(C, 9, NULL, 9, NULL);
 66: #endif

 68:   /* forms the element stiffness for the Laplacian */
 69:   FormElementStiffness(h * h, Ke);

 71:   /* assemble the matrix: only process 0 adds the values, not scalable */
 72:   if (rank == 0) {
 73:     for (i = 0; i < M; i++) {
 74:       /* node numbers for the four corners of element */
 75:       idx[0] = (m + 1) * (i / m) + (i % m);
 76:       idx[1] = idx[0] + 1;
 77:       idx[2] = idx[1] + m + 1;
 78:       idx[3] = idx[2] - 1;
 79:       if (i == M - 1 && !usezerorows) { /* If MatZeroRows not supported -> make it non-singular */
 80:         for (PetscInt ii = 0; ii < 4; ii++) {
 81:           Ke[2 * 4 + ii] = 0.0;
 82:           Ke[ii * 4 + 2] = 0.0;
 83:         }
 84:         Ke[10] = 1.0;
 85:       }
 86:       MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
 87:     }
 88:   }
 89:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
 90:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

 92:   /* create right hand side and solution */
 93:   MatCreateVecs(C, &u, &b);
 94:   VecDuplicate(u, &ustar);
 95:   VecSet(u, 0.0);
 96:   VecSet(b, 0.0);

 98:   /* assemble the right hand side: only process 0 adds the values, not scalable */
 99:   if (rank == 0) {
100:     for (i = 0; i < M; i++) {
101:       /* location of lower left corner of element */
102:       x = h * (i % m);
103:       y = h * (i / m);
104:       /* node numbers for the four corners of element */
105:       idx[0] = (m + 1) * (i / m) + (i % m);
106:       idx[1] = idx[0] + 1;
107:       idx[2] = idx[1] + m + 1;
108:       idx[3] = idx[2] - 1;
109:       FormElementRhs(x, y, h * h, r);
110:       VecSetValues(b, 4, idx, r, ADD_VALUES);
111:     }
112:   }
113:   VecAssemblyBegin(b);
114:   VecAssemblyEnd(b);

116:   /* modify matrix and rhs for Dirichlet boundary conditions */
117:   if (rank == 0) {
118:     PetscMalloc1(4 * m + 1, &rows);
119:     for (i = 0; i < m + 1; i++) {
120:       rows[i]             = i;               /* bottom */
121:       rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
122:     }
123:     count = m + 1; /* left side */
124:     for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;

126:     count = 2 * m; /* right side */
127:     for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;

129:     for (i = 0; i < 4 * m; i++) {
130:       val = h * (rows[i] / (m + 1));
131:       VecSetValues(u, 1, &rows[i], &val, INSERT_VALUES);
132:       VecSetValues(b, 1, &rows[i], &val, INSERT_VALUES);
133:     }
134:     if (usezerorows) MatZeroRows(C, 4 * m, rows, 1.0, 0, 0);
135:     PetscFree(rows);
136:   } else {
137:     if (usezerorows) MatZeroRows(C, 0, NULL, 1.0, 0, 0);
138:   }
139:   VecAssemblyBegin(u);
140:   VecAssemblyEnd(u);
141:   VecAssemblyBegin(b);
142:   VecAssemblyEnd(b);

144:   if (!usezerorows) {
145:     VecSet(ustar, 1.0);
146:     MatMult(C, ustar, b);
147:   }

149:   /* solve linear system */
150:   KSPCreate(PETSC_COMM_WORLD, &ksp);
151:   KSPSetOperators(ksp, C, C);
152:   KSPSetFromOptions(ksp);
153:   KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
154:   KSPSolve(ksp, b, u);

156:   /* check error */
157:   if (usezerorows) {
158:     if (rank == 0) {
159:       for (i = 0; i < N; i++) {
160:         val = h * (i / (m + 1));
161:         VecSetValues(ustar, 1, &i, &val, INSERT_VALUES);
162:       }
163:     }
164:     VecAssemblyBegin(ustar);
165:     VecAssemblyEnd(ustar);
166:   }

168:   VecAXPY(u, -1.0, ustar);
169:   VecNorm(u, NORM_2, &norm);
170:   KSPGetIterationNumber(ksp, &its);
171:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its);

173:   KSPDestroy(&ksp);
174:   VecDestroy(&ustar);
175:   VecDestroy(&u);
176:   VecDestroy(&b);
177:   MatDestroy(&C);
178:   PetscFinalize();
179:   return 0;
180: }

182: /*TEST

184:     test:
185:       args: -ksp_monitor_short -m 5 -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

187:     test:
188:       suffix: 3
189:       args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always

191:     test:
192:       suffix: 5
193:       args: -pc_type eisenstat -ksp_monitor_short -m 5  -ksp_gmres_cgs_refinement_type refine_always

195:     test:
196:       requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
197:       suffix: hypre_device_none
198:       output_file: output/ex4_hypre_none.out
199:       nsize: {{1 2}}
200:       args: -usezerorows 0 -mat_type hypre -pc_type none -m 5

202:     test:
203:       requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
204:       suffix: hypre_device_amg
205:       nsize: {{1 2}}
206:       args: -usezerorows 0 -mat_type hypre -pc_type hypre -m 25 -ksp_type cg -ksp_norm_type natural

208: TEST*/