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: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
 12: {
 13:   PetscFunctionBeginUser;
 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:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

 33: PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
 34: {
 35:   PetscFunctionBeginUser;
 36:   r[0] = 0.;
 37:   r[1] = 0.;
 38:   r[2] = 0.;
 39:   r[3] = 0.0;
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

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

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

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

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

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

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

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

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

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

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

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

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

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

159:   /* Check error */
160:   PetscCall(VecGetOwnershipRange(ustar, &start, &end));
161:   for (i = start; i < end; i++) {
162:     val = h * (i / (m + 1));
163:     PetscCall(VecSetValues(ustar, 1, &i, &val, INSERT_VALUES));
164:   }
165:   PetscCall(VecAssemblyBegin(ustar));
166:   PetscCall(VecAssemblyEnd(ustar));
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:   /* Free work space */
173:   PetscCall(KSPDestroy(&ksp));
174:   PetscCall(VecDestroy(&ustar));
175:   PetscCall(VecDestroy(&u));
176:   PetscCall(VecDestroy(&b));
177:   PetscCall(MatDestroy(&C));
178:   PetscCall(PetscFinalize());
179:   return 0;
180: }

182: /*TEST

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

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

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

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

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

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

212: TEST*/