Actual source code: ex4.c

  1: static char help[] = "Bilinear elements on the unit square for the Laplacian. Input arguments are:\n\
  2:   -m <size> : problem size\n\n";

  4: #include <petscksp.h>

  6: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
  7: {
  8:   Ke[0]  = H / 6.0;
  9:   Ke[1]  = -.125 * H;
 10:   Ke[2]  = H / 12.0;
 11:   Ke[3]  = -.125 * H;
 12:   Ke[4]  = -.125 * H;
 13:   Ke[5]  = H / 6.0;
 14:   Ke[6]  = -.125 * H;
 15:   Ke[7]  = H / 12.0;
 16:   Ke[8]  = H / 12.0;
 17:   Ke[9]  = -.125 * H;
 18:   Ke[10] = H / 6.0;
 19:   Ke[11] = -.125 * H;
 20:   Ke[12] = -.125 * H;
 21:   Ke[13] = H / 12.0;
 22:   Ke[14] = -.125 * H;
 23:   Ke[15] = H / 6.0;
 24:   return PETSC_SUCCESS;
 25: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

172:   { // Test getting Jacobi diag
173:     PC        pc;
174:     PetscBool is_pcjacobi;

176:     PetscCall(KSPGetPC(ksp, &pc));
177:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_pcjacobi));
178:     if (is_pcjacobi) {
179:       Vec diag;

181:       PetscCall(MatCreateVecs(C, &diag, NULL));
182:       PetscCall(PCJacobiGetDiagonal(pc, diag, NULL));
183:       PetscCall(VecNorm(diag, NORM_2, &norm));
184:       PetscCheck(norm > 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Jacobi preconditioner should have norm greater than 0");
185:       PetscCall(VecDestroy(&diag));
186:     }
187:   }

189:   PetscCall(KSPDestroy(&ksp));
190:   PetscCall(VecDestroy(&ustar));
191:   PetscCall(VecDestroy(&u));
192:   PetscCall(VecDestroy(&b));
193:   PetscCall(MatDestroy(&C));
194:   PetscCall(PetscFinalize());
195:   return 0;
196: }

198: /*TEST

200:     test:
201:       args: -ksp_monitor_short -m 5 -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always

203:     test:
204:       suffix: 3
205:       args: -pc_type sor -pc_sor_symmetric -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always

207:     test:
208:       suffix: 5
209:       args: -pc_type eisenstat -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always

211:     test:
212:       requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
213:       suffix: hypre_device_none
214:       output_file: output/ex4_hypre_none.out
215:       nsize: {{1 2}}
216:       args: -usezerorows 0 -mat_type hypre -pc_type none -m 5

218:     test:
219:       requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
220:       suffix: hypre_device_amg
221:       nsize: {{1 2}}
222:       args: -usezerorows 0 -mat_type hypre -pc_type hypre -m 25 -ksp_type cg -ksp_norm_type natural

224:     test:
225:       requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
226:       suffix: dev
227:       nsize: 2
228:       args: -m 25 -mat_type hypre -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_coarsen_type PMIS

230: TEST*/