1 static char help[] = "Test TSComputeIJacobian()\n\n"; 2 3 #include <petscsys.h> 4 #include <petscdm.h> 5 #include <petscdmda.h> 6 #include <petscts.h> 7 8 typedef struct { 9 PetscScalar u, v; 10 } Field; 11 12 typedef struct { 13 PetscReal D1, D2, gamma, kappa; 14 } AppCtx; 15 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]; 27 28 PetscFunctionBeginUser; 29 PetscCall(TSGetDM(ts, &da)); 30 PetscCall(DMGetLocalVector(da, &localU)); 31 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)); 32 33 hx = 2.50 / (PetscReal)Mx; 34 sx = 1.0 / (hx * hx); 35 hy = 2.50 / (PetscReal)My; 36 sy = 1.0 / (hy * hy); 37 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 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU)); 45 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU)); 46 47 /* 48 Get pointers to vector data 49 */ 50 PetscCall(DMDAVecGetArrayRead(da, localU, &u)); 51 52 /* 53 Get local grid boundaries 54 */ 55 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 56 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; 79 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; 82 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);*/ 86 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; 107 108 PetscCall(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; 122 123 PetscCall(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 } 127 128 /* 129 Restore vectors 130 */ 131 PetscCall(PetscLogFlops(19 * xm * ym)); 132 PetscCall(DMDAVecRestoreArrayRead(da, localU, &u)); 133 PetscCall(DMRestoreLocalVector(da, &localU)); 134 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 135 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 136 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 137 PetscFunctionReturn(0); 138 } 139 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; 145 146 PetscFunctionBeginUser; 147 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)); 148 149 hx = 2.5 / (PetscReal)Mx; 150 hy = 2.5 / (PetscReal)My; 151 152 /* 153 Get pointers to vector data 154 */ 155 PetscCall(DMDAVecGetArray(da, U, &u)); 156 157 /* 158 Get local grid boundaries 159 */ 160 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 161 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; 171 172 u[j][i].u = 1.0 - 2.0 * u[j][i].v; 173 } 174 } 175 176 /* 177 Restore vectors 178 */ 179 PetscCall(DMDAVecRestoreArray(da, U, &u)); 180 PetscFunctionReturn(0); 181 } 182 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; 191 192 PetscFunctionBeginUser; 193 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 194 195 appctx.D1 = 8.0e-5; 196 appctx.D2 = 4.0e-5; 197 appctx.gamma = .024; 198 appctx.kappa = .06; 199 200 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 201 Create distributed array (DMDA) to manage parallel grid and vectors 202 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 203 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)); 204 PetscCall(DMSetFromOptions(da)); 205 PetscCall(DMSetUp(da)); 206 PetscCall(DMDASetFieldName(da, 0, "u")); 207 PetscCall(DMDASetFieldName(da, 1, "v")); 208 PetscCall(DMSetMatType(da, MATAIJ)); 209 PetscCall(DMCreateMatrix(da, &Jac)); 210 PetscCall(DMCreateMatrix(da, &Jac2)); 211 212 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 213 Extract global vectors from DMDA; then duplicate for remaining 214 vectors that are the same types 215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 216 PetscCall(DMCreateGlobalVector(da, &U)); 217 PetscCall(VecDuplicate(U, &Udot)); 218 PetscCall(VecSet(Udot, 0.0)); 219 220 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 221 Create timestepping solver context 222 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 223 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 224 PetscCall(TSSetDM(ts, da)); 225 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 226 PetscCall(TSSetRHSJacobian(ts, Jac, Jac, RHSJacobian, &appctx)); 227 228 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 229 Set initial conditions 230 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 231 PetscCall(InitialConditions(da, U)); 232 PetscCall(TSSetSolution(ts, U)); 233 PetscCall(TSSetFromOptions(ts)); 234 PetscCall(TSSetUp(ts)); 235 236 shift = 2.; 237 PetscCall(TSComputeIJacobian(ts, t, U, Udot, shift, Jac2, Jac2, PETSC_FALSE)); 238 shift = 1.; 239 PetscCall(TSComputeIJacobian(ts, t, U, Udot, shift, Jac, Jac, PETSC_FALSE)); 240 shift = 2.; 241 PetscCall(TSComputeIJacobian(ts, t, U, Udot, shift, Jac, Jac, PETSC_FALSE)); 242 PetscCall(MatAXPY(Jac, -1, Jac2, SAME_NONZERO_PATTERN)); 243 PetscCall(MatNorm(Jac, NORM_INFINITY, &norm)); 244 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)); 245 PetscCall(MatDestroy(&Jac)); 246 PetscCall(MatDestroy(&Jac2)); 247 PetscCall(VecDestroy(&U)); 248 PetscCall(VecDestroy(&Udot)); 249 PetscCall(TSDestroy(&ts)); 250 PetscCall(DMDestroy(&da)); 251 PetscCall(PetscFinalize()); 252 return 0; 253 } 254 255 /*TEST 256 test: 257 258 TEST*/ 259