Actual source code: ex3.c

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

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

  8: #include <petscksp.h>

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

 32: PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
 33: {
 34:   PetscFunctionBeginUser;
 35:   r[0] = 0.;
 36:   r[1] = 0.;
 37:   r[2] = 0.;
 38:   r[3] = 0.0;
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 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, build_sol;
 51:   KSP         ksp;
 52:   PetscBool   viewkspest = PETSC_FALSE, testbuildsolution = PETSC_FALSE;

 54:   PetscFunctionBeginUser;
 55:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 56:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 57:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_est_view", &viewkspest, NULL));
 58:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_build_solution", &testbuildsolution, 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 (testbuildsolution) {
153:     PetscBool ok;

155:     PetscCall(VecDuplicate(u, &build_sol));
156:     PetscCall(KSPBuildSolution(ksp, build_sol, NULL));
157:     PetscCall(VecEqual(u, build_sol, &ok));
158:     PetscCheck(ok, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "KSPBuildSolution() returned incorrect solution");
159:     PetscCall(VecDestroy(&build_sol));
160:   }

162:   if (viewkspest) {
163:     KSP kspest;

165:     PetscCall(KSPChebyshevEstEigGetKSP(ksp, &kspest));
166:     if (kspest) PetscCall(KSPView(kspest, PETSC_VIEWER_STDOUT_WORLD));
167:   }

169:   /* Check error */
170:   PetscCall(VecGetOwnershipRange(ustar, &start, &end));
171:   for (i = start; i < end; i++) {
172:     val = h * (i / (m + 1));
173:     PetscCall(VecSetValues(ustar, 1, &i, &val, INSERT_VALUES));
174:   }
175:   PetscCall(VecAssemblyBegin(ustar));
176:   PetscCall(VecAssemblyEnd(ustar));
177:   PetscCall(VecAXPY(u, -1.0, ustar));
178:   PetscCall(VecNorm(u, NORM_2, &norm));
179:   PetscCall(KSPGetIterationNumber(ksp, &its));
180:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its));

182:   /* Free work space */
183:   PetscCall(KSPDestroy(&ksp));
184:   PetscCall(VecDestroy(&ustar));
185:   PetscCall(VecDestroy(&u));
186:   PetscCall(VecDestroy(&b));
187:   PetscCall(MatDestroy(&C));
188:   PetscCall(PetscFinalize());
189:   return 0;
190: }

192: /*TEST

194:     test:
195:       args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always

197:     test:
198:       suffix: 2
199:       nsize: 2
200:       args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always

202:     test:
203:       suffix: 2_kokkos
204:       nsize: 2
205:       args: -vec_mdot_use_gemv {{0 1}} -vec_maxpy_use_gemv {{0 1}}
206:       args: -pc_type jacobi -ksp_monitor_short -m 5 -ksp_gmres_cgs_refinement_type refine_always -mat_type aijkokkos -vec_type kokkos
207:       output_file: output/ex3_2.out
208:       requires: kokkos_kernels

210:     test:
211:       suffix: nocheby
212:       args: -ksp_est_view

214:     test:
215:       suffix: chebynoest
216:       args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_eigenvalues 0.1,1.0

218:     test:
219:       suffix: chebyest
220:       args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_esteig
221:       filter:  sed -e "s/Iterations 19/Iterations 20/g"

223:     test:
224:       suffix: gamg_provided_not_ok
225:       filter: grep -v "variant HERMITIAN" | sed -e "s/Iterations 4/Iterations 5/g"
226:       args: -pc_type gamg -mg_levels_pc_type sor -mg_levels_esteig_ksp_type cg -ksp_view

228:     test:
229:       suffix: build_solution
230:       requires: !complex
231:       filter: grep -v Norm
232:       args: -ksp_type {{chebyshev cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr bicg minres lcd gcr cgls richardson}} -test_build_solution

234: TEST*/