Actual source code: ex24.c

  1: static char help[] = "Test TSComputeIJacobian()\n\n";

  3: #include <petscsys.h>
  4: #include <petscdm.h>
  5: #include <petscdmda.h>
  6: #include <petscts.h>

  8: typedef struct {
  9:   PetscScalar u, v;
 10: } Field;

 12: typedef struct {
 13:   PetscReal D1, D2, gamma, kappa;
 14: } AppCtx;

 16: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat BB, void *ctx)
 17: {
 18:   AppCtx     *appctx = (AppCtx *)ctx; /* user-defined application context */
 19:   DM          da;
 20:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
 21:   PetscReal   hx, hy, sx, sy;
 22:   PetscScalar uc, vc;
 23:   Field     **u;
 24:   Vec         localU;
 25:   MatStencil  stencil[6], rowstencil;
 26:   PetscScalar entries[6];

 29:   TSGetDM(ts, &da);
 30:   DMGetLocalVector(da, &localU);
 31:   DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);

 33:   hx = 2.50 / (PetscReal)Mx;
 34:   sx = 1.0 / (hx * hx);
 35:   hy = 2.50 / (PetscReal)My;
 36:   sy = 1.0 / (hy * hy);

 38:   /*
 39:      Scatter ghost points to local vector,using the 2-step process
 40:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
 41:      By placing code between these two statements, computations can be
 42:      done while messages are in transition.
 43:   */
 44:   DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
 45:   DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);

 47:   /*
 48:      Get pointers to vector data
 49:   */
 50:   DMDAVecGetArrayRead(da, localU, &u);

 52:   /*
 53:      Get local grid boundaries
 54:   */
 55:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

 57:   stencil[0].k = 0;
 58:   stencil[1].k = 0;
 59:   stencil[2].k = 0;
 60:   stencil[3].k = 0;
 61:   stencil[4].k = 0;
 62:   stencil[5].k = 0;
 63:   rowstencil.k = 0;
 64:   /*
 65:      Compute function over the locally owned part of the grid
 66:   */
 67:   for (j = ys; j < ys + ym; j++) {
 68:     stencil[0].j = j - 1;
 69:     stencil[1].j = j + 1;
 70:     stencil[2].j = j;
 71:     stencil[3].j = j;
 72:     stencil[4].j = j;
 73:     stencil[5].j = j;
 74:     rowstencil.k = 0;
 75:     rowstencil.j = j;
 76:     for (i = xs; i < xs + xm; i++) {
 77:       uc = u[j][i].u;
 78:       vc = u[j][i].v;

 80:       /*      uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
 81:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;

 83:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
 84:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
 85:        f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/

 87:       stencil[0].i = i;
 88:       stencil[0].c = 0;
 89:       entries[0]   = appctx->D1 * sy;
 90:       stencil[1].i = i;
 91:       stencil[1].c = 0;
 92:       entries[1]   = appctx->D1 * sy;
 93:       stencil[2].i = i - 1;
 94:       stencil[2].c = 0;
 95:       entries[2]   = appctx->D1 * sx;
 96:       stencil[3].i = i + 1;
 97:       stencil[3].c = 0;
 98:       entries[3]   = appctx->D1 * sx;
 99:       stencil[4].i = i;
100:       stencil[4].c = 0;
101:       entries[4]   = -2.0 * appctx->D1 * (sx + sy) - vc * vc - appctx->gamma;
102:       stencil[5].i = i;
103:       stencil[5].c = 1;
104:       entries[5]   = -2.0 * uc * vc;
105:       rowstencil.i = i;
106:       rowstencil.c = 0;

108:       MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
109:       stencil[0].c = 1;
110:       entries[0]   = appctx->D2 * sy;
111:       stencil[1].c = 1;
112:       entries[1]   = appctx->D2 * sy;
113:       stencil[2].c = 1;
114:       entries[2]   = appctx->D2 * sx;
115:       stencil[3].c = 1;
116:       entries[3]   = appctx->D2 * sx;
117:       stencil[4].c = 1;
118:       entries[4]   = -2.0 * appctx->D2 * (sx + sy) + 2.0 * uc * vc - appctx->gamma - appctx->kappa;
119:       stencil[5].c = 0;
120:       entries[5]   = vc * vc;
121:       rowstencil.c = 1;

123:       MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES);
124:       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
125:     }
126:   }

128:   /*
129:      Restore vectors
130:   */
131:   PetscLogFlops(19 * xm * ym);
132:   DMDAVecRestoreArrayRead(da, localU, &u);
133:   DMRestoreLocalVector(da, &localU);
134:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
135:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
136:   MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
137:   return 0;
138: }

140: PetscErrorCode InitialConditions(DM da, Vec U)
141: {
142:   PetscInt  i, j, xs, ys, xm, ym, Mx, My;
143:   Field   **u;
144:   PetscReal hx, hy, x, y;

147:   DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);

149:   hx = 2.5 / (PetscReal)Mx;
150:   hy = 2.5 / (PetscReal)My;

152:   /*
153:      Get pointers to vector data
154:   */
155:   DMDAVecGetArray(da, U, &u);

157:   /*
158:      Get local grid boundaries
159:   */
160:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

162:   /*
163:      Compute function over the locally owned part of the grid
164:   */
165:   for (j = ys; j < ys + ym; j++) {
166:     y = j * hy;
167:     for (i = xs; i < xs + xm; i++) {
168:       x = i * hx;
169:       if ((1.0 <= x) && (x <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25 * PetscPowReal(PetscSinReal(4.0 * PETSC_PI * x), 2.0) * PetscPowReal(PetscSinReal(4.0 * PETSC_PI * y), 2.0);
170:       else u[j][i].v = 0.0;

172:       u[j][i].u = 1.0 - 2.0 * u[j][i].v;
173:     }
174:   }

176:   /*
177:      Restore vectors
178:   */
179:   DMDAVecRestoreArray(da, U, &u);
180:   return 0;
181: }

183: int main(int argc, char **argv)
184: {
185:   TS        ts;
186:   Vec       U, Udot;
187:   Mat       Jac, Jac2;
188:   DM        da;
189:   AppCtx    appctx;
190:   PetscReal t = 0, shift, norm;

193:   PetscInitialize(&argc, &argv, (char *)0, help);

195:   appctx.D1    = 8.0e-5;
196:   appctx.D2    = 4.0e-5;
197:   appctx.gamma = .024;
198:   appctx.kappa = .06;

200:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201:      Create distributed array (DMDA) to manage parallel grid and vectors
202:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 64, 64, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, NULL, &da);
204:   DMSetFromOptions(da);
205:   DMSetUp(da);
206:   DMDASetFieldName(da, 0, "u");
207:   DMDASetFieldName(da, 1, "v");
208:   DMSetMatType(da, MATAIJ);
209:   DMCreateMatrix(da, &Jac);
210:   DMCreateMatrix(da, &Jac2);

212:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213:      Extract global vectors from DMDA; then duplicate for remaining
214:      vectors that are the same types
215:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
216:   DMCreateGlobalVector(da, &U);
217:   VecDuplicate(U, &Udot);
218:   VecSet(Udot, 0.0);

220:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221:      Create timestepping solver context
222:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
223:   TSCreate(PETSC_COMM_WORLD, &ts);
224:   TSSetDM(ts, da);
225:   TSSetProblemType(ts, TS_NONLINEAR);
226:   TSSetRHSJacobian(ts, Jac, Jac, RHSJacobian, &appctx);

228:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229:      Set initial conditions
230:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
231:   InitialConditions(da, U);
232:   TSSetSolution(ts, U);
233:   TSSetFromOptions(ts);
234:   TSSetUp(ts);

236:   shift = 2.;
237:   TSComputeIJacobian(ts, t, U, Udot, shift, Jac2, Jac2, PETSC_FALSE);
238:   shift = 1.;
239:   TSComputeIJacobian(ts, t, U, Udot, shift, Jac, Jac, PETSC_FALSE);
240:   shift = 2.;
241:   TSComputeIJacobian(ts, t, U, Udot, shift, Jac, Jac, PETSC_FALSE);
242:   MatAXPY(Jac, -1, Jac2, SAME_NONZERO_PATTERN);
243:   MatNorm(Jac, NORM_INFINITY, &norm);
244:   if (norm > 100.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Error Norm %g \n Incorrect behaviour of TSComputeIJacobian(). The two matrices should have the same results.\n", (double)norm);
245:   MatDestroy(&Jac);
246:   MatDestroy(&Jac2);
247:   VecDestroy(&U);
248:   VecDestroy(&Udot);
249:   TSDestroy(&ts);
250:   DMDestroy(&da);
251:   PetscFinalize();
252:   return 0;
253: }

255: /*TEST
256:   test:

258: TEST*/