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;

 31:   PetscInitialize(&argc, &argv, (char *)0, help);
 32:   PetscOptionsGetInt(NULL, NULL, "-dof", &dof, NULL);
 33:   PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL);
 34:   PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL);

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

 48:   DMCreateGlobalVector(da, &x);
 49:   DMCreateGlobalVector(da, &b);
 50:   ComputeRHS(da, b);
 51:   DMSetMatType(da, MATBAIJ);
 52:   DMSetFromOptions(da);
 53:   DMCreateMatrix(da, &A);
 54:   ComputeMatrix(da, A);

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

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

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

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

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

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

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

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

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

132:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
133:   /* For simplicity, this example only works on mx=my=mz */

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

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

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

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

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

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

215: /*TEST

217:    test:
218:       args: -ksp_monitor_short -dm_mat_type sbaij -ksp_monitor_short -pc_type cholesky -ksp_view

220: TEST*/