Actual source code: ex3.c


  2: static char help[] = "Bilinear elements on the unit square for Laplacian.  To test the parallel\n\
  3: matrix assembly, the matrix is intentionally laid out across processors\n\
  4: differently from the way it is assembled.  Input arguments are:\n\
  5:   -m <size> : problem size\n\n";

  7: /* Addendum: piggy-backing on this example to test KSPChebyshev methods */

  9: #include <petscksp.h>

 11: int FormElementStiffness(PetscReal H, PetscScalar *Ke)
 12: {
 14:   Ke[0]  = H / 6.0;
 15:   Ke[1]  = -.125 * H;
 16:   Ke[2]  = H / 12.0;
 17:   Ke[3]  = -.125 * H;
 18:   Ke[4]  = -.125 * H;
 19:   Ke[5]  = H / 6.0;
 20:   Ke[6]  = -.125 * H;
 21:   Ke[7]  = H / 12.0;
 22:   Ke[8]  = H / 12.0;
 23:   Ke[9]  = -.125 * H;
 24:   Ke[10] = H / 6.0;
 25:   Ke[11] = -.125 * H;
 26:   Ke[12] = -.125 * H;
 27:   Ke[13] = H / 12.0;
 28:   Ke[14] = -.125 * H;
 29:   Ke[15] = H / 6.0;
 30:   return 0;
 31: }
 32: int FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
 33: {
 35:   r[0] = 0.;
 36:   r[1] = 0.;
 37:   r[2] = 0.;
 38:   r[3] = 0.0;
 39:   return 0;
 40: }

 42: int main(int argc, char **args)
 43: {
 44:   Mat         C;
 45:   PetscMPIInt rank, size;
 46:   PetscInt    i, m = 5, N, start, end, M, its;
 47:   PetscScalar val, Ke[16], r[4];
 48:   PetscReal   x, y, h, norm;
 49:   PetscInt    idx[4], count, *rows;
 50:   Vec         u, ustar, b;
 51:   KSP         ksp;
 52:   PetscBool   viewkspest = PETSC_FALSE;

 55:   PetscInitialize(&argc, &args, (char *)0, help);
 56:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
 57:   PetscOptionsGetBool(NULL, NULL, "-ksp_est_view", &viewkspest, NULL);
 58:   N = (m + 1) * (m + 1); /* dimension of matrix */
 59:   M = m * m;             /* number of elements */
 60:   h = 1.0 / m;           /* mesh width */
 61:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 62:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 64:   /* Create stiffness matrix */
 65:   MatCreate(PETSC_COMM_WORLD, &C);
 66:   MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N);
 67:   MatSetFromOptions(C);
 68:   MatSetUp(C);
 69:   start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
 70:   end   = start + M / size + ((M % size) > rank);

 72:   /* Assemble matrix */
 73:   FormElementStiffness(h * h, Ke); /* element stiffness for Laplacian */
 74:   for (i = start; i < end; i++) {
 75:     /* node numbers for the four corners of element */
 76:     idx[0] = (m + 1) * (i / m) + (i % m);
 77:     idx[1] = idx[0] + 1;
 78:     idx[2] = idx[1] + m + 1;
 79:     idx[3] = idx[2] - 1;
 80:     MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
 81:   }
 82:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

 85:   /* Create right-hand-side and solution vectors */
 86:   VecCreate(PETSC_COMM_WORLD, &u);
 87:   VecSetSizes(u, PETSC_DECIDE, N);
 88:   VecSetFromOptions(u);
 89:   PetscObjectSetName((PetscObject)u, "Approx. Solution");
 90:   VecDuplicate(u, &b);
 91:   PetscObjectSetName((PetscObject)b, "Right hand side");
 92:   VecDuplicate(b, &ustar);
 93:   VecSet(u, 0.0);
 94:   VecSet(b, 0.0);

 96:   /* Assemble right-hand-side vector */
 97:   for (i = start; i < end; i++) {
 98:     /* location of lower left corner of element */
 99:     x = h * (i % m);
100:     y = h * (i / m);
101:     /* node numbers for the four corners of element */
102:     idx[0] = (m + 1) * (i / m) + (i % m);
103:     idx[1] = idx[0] + 1;
104:     idx[2] = idx[1] + m + 1;
105:     idx[3] = idx[2] - 1;
106:     FormElementRhs(x, y, h * h, r);
107:     VecSetValues(b, 4, idx, r, ADD_VALUES);
108:   }
109:   VecAssemblyBegin(b);
110:   VecAssemblyEnd(b);

112:   /* Modify matrix and right-hand-side for Dirichlet boundary conditions */
113:   PetscMalloc1(4 * m, &rows);
114:   for (i = 0; i < m + 1; i++) {
115:     rows[i]             = i;               /* bottom */
116:     rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
117:   }
118:   count = m + 1; /* left side */
119:   for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;

121:   count = 2 * m; /* left side */
122:   for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
123:   for (i = 0; i < 4 * m; i++) {
124:     val = h * (rows[i] / (m + 1));
125:     VecSetValues(u, 1, &rows[i], &val, INSERT_VALUES);
126:     VecSetValues(b, 1, &rows[i], &val, INSERT_VALUES);
127:   }
128:   MatZeroRows(C, 4 * m, rows, 1.0, 0, 0);

130:   PetscFree(rows);
131:   VecAssemblyBegin(u);
132:   VecAssemblyEnd(u);
133:   VecAssemblyBegin(b);
134:   VecAssemblyEnd(b);

136:   {
137:     Mat A;
138:     MatConvert(C, MATSAME, MAT_INITIAL_MATRIX, &A);
139:     MatDestroy(&C);
140:     MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C);
141:     MatDestroy(&A);
142:   }

144:   /* Solve linear system */
145:   KSPCreate(PETSC_COMM_WORLD, &ksp);
146:   KSPSetOperators(ksp, C, C);
147:   KSPSetFromOptions(ksp);
148:   KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
149:   KSPSolve(ksp, b, u);

151:   if (viewkspest) {
152:     KSP kspest;

154:     KSPChebyshevEstEigGetKSP(ksp, &kspest);
155:     if (kspest) KSPView(kspest, PETSC_VIEWER_STDOUT_WORLD);
156:   }

158:   /* Check error */
159:   VecGetOwnershipRange(ustar, &start, &end);
160:   for (i = start; i < end; i++) {
161:     val = h * (i / (m + 1));
162:     VecSetValues(ustar, 1, &i, &val, INSERT_VALUES);
163:   }
164:   VecAssemblyBegin(ustar);
165:   VecAssemblyEnd(ustar);
166:   VecAXPY(u, -1.0, ustar);
167:   VecNorm(u, NORM_2, &norm);
168:   KSPGetIterationNumber(ksp, &its);
169:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its);

171:   /* Free work space */
172:   KSPDestroy(&ksp);
173:   VecDestroy(&ustar);
174:   VecDestroy(&u);
175:   VecDestroy(&b);
176:   MatDestroy(&C);
177:   PetscFinalize();
178:   return 0;
179: }

181: /*TEST

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

186:     test:
187:       suffix: 2
188:       nsize: 2
189:       args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always

191:     test:
192:       suffix: 2_kokkos
193:       nsize: 2
194:       args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always -mat_type aijkokkos -vec_type kokkos
195:       output_file: output/ex3_2.out
196:       requires: kokkos_kernels

198:     test:
199:       suffix: nocheby
200:       args: -ksp_est_view

202:     test:
203:       suffix: chebynoest
204:       args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_eigenvalues 0.1,1.0

206:     test:
207:       suffix: chebyest
208:       args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_esteig
209:       filter:  sed -e "s/Iterations 19/Iterations 20/g"

211: TEST*/