Actual source code: ex32.c

  1: /*
  2:   Laplacian in 3D. Use for testing BAIJ matrix.
  3:   Modeled by the partial differential equation

  5:    - Laplacian u = 1,0 < x,y,z < 1,

  7:    with boundary conditions
  8:    u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
  9: */

 11: static char help[] = "Solves 3D Laplacian using wirebasket based multigrid.\n\n";

 13: #include <petscdm.h>
 14: #include <petscdmda.h>
 15: #include <petscksp.h>

 17: extern PetscErrorCode ComputeMatrix(DM, Mat);
 18: extern PetscErrorCode ComputeRHS(DM, Vec);

 20: int main(int argc, char **argv)
 21: {
 22:   KSP       ksp;
 23:   PC        pc;
 24:   Vec       x, b;
 25:   DM        da;
 26:   Mat       A, Atrans;
 27:   PetscInt  dof = 1, M = 8;
 28:   PetscBool flg, trans = PETSC_FALSE, sbaij = PETSC_FALSE;

 30:   PetscFunctionBeginUser;
 31:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 32:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL));
 33:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
 34:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL));
 35:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-sbaij", &sbaij, NULL));

 37:   PetscCall(DMDACreate(PETSC_COMM_WORLD, &da));
 38:   PetscCall(DMSetDimension(da, 3));
 39:   PetscCall(DMDASetBoundaryType(da, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE));
 40:   PetscCall(DMDASetStencilType(da, DMDA_STENCIL_STAR));
 41:   PetscCall(DMDASetSizes(da, M, M, M));
 42:   PetscCall(DMDASetNumProcs(da, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
 43:   PetscCall(DMDASetDof(da, dof));
 44:   PetscCall(DMDASetStencilWidth(da, 1));
 45:   PetscCall(DMDASetOwnershipRanges(da, NULL, NULL, NULL));
 46:   PetscCall(DMSetFromOptions(da));
 47:   PetscCall(DMSetUp(da));

 49:   PetscCall(DMCreateGlobalVector(da, &x));
 50:   PetscCall(DMCreateGlobalVector(da, &b));
 51:   PetscCall(ComputeRHS(da, b));
 52:   PetscCall(DMSetMatType(da, MATBAIJ));
 53:   if (sbaij) PetscCall(DMSetMatType(da, MATSBAIJ));
 54:   PetscCall(DMCreateMatrix(da, &A));
 55:   PetscCall(ComputeMatrix(da, A));

 57:   /* A is non-symmetric. Make A = 0.5*(A + Atrans) symmetric for testing icc and cholesky */
 58:   PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
 59:   PetscCall(MatAXPY(A, 1.0, Atrans, DIFFERENT_NONZERO_PATTERN));
 60:   PetscCall(MatScale(A, 0.5));
 61:   PetscCall(MatDestroy(&Atrans));

 63:   /* Test sbaij matrix */
 64:   flg = PETSC_FALSE;
 65:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_sbaij1", &flg, NULL));
 66:   if (flg) {
 67:     Mat       sA;
 68:     PetscBool issymm;
 69:     PetscCall(MatIsTranspose(A, A, 0.0, &issymm));
 70:     if (issymm) {
 71:       PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
 72:     } else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: A is non-symmetric\n"));
 73:     PetscCall(MatConvert(A, MATSBAIJ, MAT_INITIAL_MATRIX, &sA));
 74:     PetscCall(MatDestroy(&A));
 75:     A = sA;
 76:   }

 78:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
 79:   PetscCall(KSPSetFromOptions(ksp));
 80:   PetscCall(KSPSetOperators(ksp, A, A));
 81:   PetscCall(KSPGetPC(ksp, &pc));
 82:   PetscCall(PCSetDM(pc, (DM)da));

 84:   if (trans) {
 85:     PetscCall(KSPSolveTranspose(ksp, b, x));
 86:   } else {
 87:     PetscCall(KSPSolve(ksp, b, x));
 88:   }

 90:   /* check final residual */
 91:   flg = PETSC_FALSE;
 92:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-check_final_residual", &flg, NULL));
 93:   if (flg) {
 94:     Vec       b1;
 95:     PetscReal norm;
 96:     PetscCall(KSPGetSolution(ksp, &x));
 97:     PetscCall(VecDuplicate(b, &b1));
 98:     PetscCall(MatMult(A, x, b1));
 99:     PetscCall(VecAXPY(b1, -1.0, b));
100:     PetscCall(VecNorm(b1, NORM_2, &norm));
101:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final residual %g\n", (double)norm));
102:     PetscCall(VecDestroy(&b1));
103:   }

105:   PetscCall(KSPDestroy(&ksp));
106:   PetscCall(VecDestroy(&x));
107:   PetscCall(VecDestroy(&b));
108:   PetscCall(MatDestroy(&A));
109:   PetscCall(DMDestroy(&da));
110:   PetscCall(PetscFinalize());
111:   return 0;
112: }

114: PetscErrorCode ComputeRHS(DM da, Vec b)
115: {
116:   PetscInt    mx, my, mz;
117:   PetscScalar h;

119:   PetscFunctionBeginUser;
120:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
121:   h = 1.0 / ((mx - 1) * (my - 1) * (mz - 1));
122:   PetscCall(VecSet(b, h));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: PetscErrorCode ComputeMatrix(DM da, Mat B)
127: {
128:   PetscInt     i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs, dof, k1, k2, k3;
129:   PetscScalar *v, *v_neighbor, Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy;
130:   MatStencil   row, col;

132:   PetscFunctionBeginUser;
133:   PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
134:   /* For simplicity, this example only works on mx=my=mz */
135:   PetscCheck(mx == my && mx == mz, PETSC_COMM_SELF, PETSC_ERR_SUP, "This example only works with mx %" PetscInt_FMT " = my %" PetscInt_FMT " = mz %" PetscInt_FMT, mx, my, mz);

137:   Hx      = 1.0 / (PetscReal)(mx - 1);
138:   Hy      = 1.0 / (PetscReal)(my - 1);
139:   Hz      = 1.0 / (PetscReal)(mz - 1);
140:   HxHydHz = Hx * Hy / Hz;
141:   HxHzdHy = Hx * Hz / Hy;
142:   HyHzdHx = Hy * Hz / Hx;

144:   PetscCall(PetscMalloc1(2 * dof * dof + 1, &v));
145:   v_neighbor = v + dof * dof;
146:   PetscCall(PetscArrayzero(v, 2 * dof * dof + 1));
147:   k3 = 0;
148:   for (k1 = 0; k1 < dof; k1++) {
149:     for (k2 = 0; k2 < dof; k2++) {
150:       if (k1 == k2) {
151:         v[k3]          = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx);
152:         v_neighbor[k3] = -HxHydHz;
153:       } else {
154:         v[k3]          = k1 / (dof * dof);
155:         v_neighbor[k3] = k2 / (dof * dof);
156:       }
157:       k3++;
158:     }
159:   }
160:   PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm));

162:   for (k = zs; k < zs + zm; k++) {
163:     for (j = ys; j < ys + ym; j++) {
164:       for (i = xs; i < xs + xm; i++) {
165:         row.i = i;
166:         row.j = j;
167:         row.k = k;
168:         if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) { /* boundary points */
169:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &row, v, INSERT_VALUES));
170:         } else { /* interior points */
171:           /* center */
172:           col.i = i;
173:           col.j = j;
174:           col.k = k;
175:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v, INSERT_VALUES));

177:           /* x neighbors */
178:           col.i = i - 1;
179:           col.j = j;
180:           col.k = k;
181:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
182:           col.i = i + 1;
183:           col.j = j;
184:           col.k = k;
185:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));

187:           /* y neighbors */
188:           col.i = i;
189:           col.j = j - 1;
190:           col.k = k;
191:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
192:           col.i = i;
193:           col.j = j + 1;
194:           col.k = k;
195:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));

197:           /* z neighbors */
198:           col.i = i;
199:           col.j = j;
200:           col.k = k - 1;
201:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
202:           col.i = i;
203:           col.j = j;
204:           col.k = k + 1;
205:           PetscCall(MatSetValuesBlockedStencil(B, 1, &row, 1, &col, v_neighbor, INSERT_VALUES));
206:         }
207:       }
208:     }
209:   }
210:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
211:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
212:   PetscCall(PetscFree(v));
213:   PetscFunctionReturn(PETSC_SUCCESS);
214: }

216: /*TEST

218:    test:
219:       args: -ksp_monitor_short -sbaij -ksp_monitor_short -pc_type cholesky -ksp_view

221: TEST*/