static char help[] = "Test TSComputeIJacobian()\n\n"; #include #include #include #include typedef struct { PetscScalar u, v; } Field; typedef struct { PetscReal D1, D2, gamma, kappa; } AppCtx; PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat BB, PetscCtx ctx) { AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */ DM da; PetscInt i, j, Mx, My, xs, ys, xm, ym; PetscReal hx, hy, sx, sy; PetscScalar uc, vc; Field **u; Vec localU; MatStencil stencil[6], rowstencil; PetscScalar entries[6]; PetscFunctionBeginUser; PetscCall(TSGetDM(ts, &da)); PetscCall(DMGetLocalVector(da, &localU)); PetscCall(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)); hx = 2.50 / (PetscReal)Mx; sx = 1.0 / (hx * hx); hy = 2.50 / (PetscReal)My; sy = 1.0 / (hy * hy); /* Scatter ghost points to local vector,using the 2-step process DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). By placing code between these two statements, computations can be done while messages are in transition. */ PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU)); PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU)); /* Get pointers to vector data */ PetscCall(DMDAVecGetArrayRead(da, localU, &u)); /* Get local grid boundaries */ PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); stencil[0].k = 0; stencil[1].k = 0; stencil[2].k = 0; stencil[3].k = 0; stencil[4].k = 0; stencil[5].k = 0; rowstencil.k = 0; /* Compute function over the locally owned part of the grid */ for (j = ys; j < ys + ym; j++) { stencil[0].j = j - 1; stencil[1].j = j + 1; stencil[2].j = j; stencil[3].j = j; stencil[4].j = j; stencil[5].j = j; rowstencil.k = 0; rowstencil.j = j; for (i = xs; i < xs + xm; i++) { uc = u[j][i].u; vc = u[j][i].v; /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/ stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1 * sy; stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1 * sy; stencil[2].i = i - 1; stencil[2].c = 0; entries[2] = appctx->D1 * sx; stencil[3].i = i + 1; stencil[3].c = 0; entries[3] = appctx->D1 * sx; stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0 * appctx->D1 * (sx + sy) - vc * vc - appctx->gamma; stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0 * uc * vc; rowstencil.i = i; rowstencil.c = 0; PetscCall(MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES)); stencil[0].c = 1; entries[0] = appctx->D2 * sy; stencil[1].c = 1; entries[1] = appctx->D2 * sy; stencil[2].c = 1; entries[2] = appctx->D2 * sx; stencil[3].c = 1; entries[3] = appctx->D2 * sx; stencil[4].c = 1; entries[4] = -2.0 * appctx->D2 * (sx + sy) + 2.0 * uc * vc - appctx->gamma - appctx->kappa; stencil[5].c = 0; entries[5] = vc * vc; rowstencil.c = 1; PetscCall(MatSetValuesStencil(A, 1, &rowstencil, 6, stencil, entries, INSERT_VALUES)); /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ } } /* Restore vectors */ PetscCall(PetscLogFlops(19 * xm * ym)); PetscCall(DMDAVecRestoreArrayRead(da, localU, &u)); PetscCall(DMRestoreLocalVector(da, &localU)); PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode InitialConditions(DM da, Vec U) { PetscInt i, j, xs, ys, xm, ym, Mx, My; Field **u; PetscReal hx, hy, x, y; PetscFunctionBeginUser; PetscCall(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)); hx = 2.5 / (PetscReal)Mx; hy = 2.5 / (PetscReal)My; /* Get pointers to vector data */ PetscCall(DMDAVecGetArray(da, U, &u)); /* Get local grid boundaries */ PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); /* Compute function over the locally owned part of the grid */ for (j = ys; j < ys + ym; j++) { y = j * hy; for (i = xs; i < xs + xm; i++) { x = i * hx; 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); else u[j][i].v = 0.0; u[j][i].u = 1.0 - 2.0 * u[j][i].v; } } /* Restore vectors */ PetscCall(DMDAVecRestoreArray(da, U, &u)); PetscFunctionReturn(PETSC_SUCCESS); } int main(int argc, char **argv) { TS ts; Vec U, Udot; Mat Jac, Jac2; DM da; AppCtx appctx; PetscReal t = 0, shift, norm; PetscFunctionBeginUser; PetscCall(PetscInitialize(&argc, &argv, NULL, help)); appctx.D1 = 8.0e-5; appctx.D2 = 4.0e-5; appctx.gamma = .024; appctx.kappa = .06; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create distributed array (DMDA) to manage parallel grid and vectors - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 64, 64, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, NULL, &da)); PetscCall(DMSetFromOptions(da)); PetscCall(DMSetUp(da)); PetscCall(DMDASetFieldName(da, 0, "u")); PetscCall(DMDASetFieldName(da, 1, "v")); PetscCall(DMSetMatType(da, MATAIJ)); PetscCall(DMCreateMatrix(da, &Jac)); PetscCall(DMCreateMatrix(da, &Jac2)); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Extract global vectors from DMDA; then duplicate for remaining vectors that are the same types - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(DMCreateGlobalVector(da, &U)); PetscCall(VecDuplicate(U, &Udot)); PetscCall(VecSet(Udot, 0.0)); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Create timestepping solver context - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); PetscCall(TSSetDM(ts, da)); PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); PetscCall(TSSetRHSJacobian(ts, Jac, Jac, RHSJacobian, &appctx)); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Set initial conditions - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(InitialConditions(da, U)); PetscCall(TSSetSolution(ts, U)); PetscCall(TSSetFromOptions(ts)); PetscCall(TSSetUp(ts)); shift = 2.; PetscCall(TSComputeIJacobian(ts, t, U, Udot, shift, Jac2, Jac2, PETSC_FALSE)); shift = 1.; PetscCall(TSComputeIJacobian(ts, t, U, Udot, shift, Jac, Jac, PETSC_FALSE)); shift = 2.; PetscCall(TSComputeIJacobian(ts, t, U, Udot, shift, Jac, Jac, PETSC_FALSE)); PetscCall(MatAXPY(Jac, -1, Jac2, SAME_NONZERO_PATTERN)); PetscCall(MatNorm(Jac, NORM_INFINITY, &norm)); if (norm > 100.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error Norm %g \n Incorrect behaviour of TSComputeIJacobian(). The two matrices should have the same results.\n", (double)norm)); PetscCall(MatDestroy(&Jac)); PetscCall(MatDestroy(&Jac2)); PetscCall(VecDestroy(&U)); PetscCall(VecDestroy(&Udot)); PetscCall(TSDestroy(&ts)); PetscCall(DMDestroy(&da)); PetscCall(PetscFinalize()); return 0; } /*TEST test: output_file: output/empty.out TEST*/