1 2 static char help[] = "Time-dependent PDE in 2d. Simplified from ex7.c for illustrating how to use TS on a structured domain. \n"; 3 /* 4 u_t = uxx + uyy 5 0 < x < 1, 0 < y < 1; 6 At t=0: u(x,y) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)) < .125 7 u(x,y) = 0.0 if r >= .125 8 9 mpiexec -n 2 ./ex13 -da_grid_x 40 -da_grid_y 40 -ts_max_steps 2 -snes_monitor -ksp_monitor 10 mpiexec -n 1 ./ex13 -snes_fd_color -ts_monitor_draw_solution 11 mpiexec -n 2 ./ex13 -ts_type sundials -ts_monitor 12 */ 13 14 #include <petscdm.h> 15 #include <petscdmda.h> 16 #include <petscts.h> 17 18 /* 19 User-defined data structures and routines 20 */ 21 typedef struct { 22 PetscReal c; 23 } AppCtx; 24 25 extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *); 26 extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *); 27 extern PetscErrorCode FormInitialSolution(DM, Vec, void *); 28 29 int main(int argc, char **argv) 30 { 31 TS ts; /* nonlinear solver */ 32 Vec u, r; /* solution, residual vector */ 33 Mat J; /* Jacobian matrix */ 34 PetscInt steps; /* iterations for convergence */ 35 DM da; 36 PetscReal ftime, dt; 37 AppCtx user; /* user-defined work context */ 38 39 PetscFunctionBeginUser; 40 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 41 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 42 Create distributed array (DMDA) to manage parallel grid and vectors 43 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 44 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 8, 8, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da)); 45 PetscCall(DMSetFromOptions(da)); 46 PetscCall(DMSetUp(da)); 47 48 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 49 Extract global vectors from DMDA; 50 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 51 PetscCall(DMCreateGlobalVector(da, &u)); 52 PetscCall(VecDuplicate(u, &r)); 53 54 /* Initialize user application context */ 55 user.c = -30.0; 56 57 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 58 Create timestepping solver context 59 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 60 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 61 PetscCall(TSSetDM(ts, da)); 62 PetscCall(TSSetType(ts, TSBEULER)); 63 PetscCall(TSSetRHSFunction(ts, r, RHSFunction, &user)); 64 65 /* Set Jacobian */ 66 PetscCall(DMSetMatType(da, MATAIJ)); 67 PetscCall(DMCreateMatrix(da, &J)); 68 PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, NULL)); 69 70 ftime = 1.0; 71 PetscCall(TSSetMaxTime(ts, ftime)); 72 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 73 74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 75 Set initial conditions 76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 77 PetscCall(FormInitialSolution(da, u, &user)); 78 dt = .01; 79 PetscCall(TSSetTimeStep(ts, dt)); 80 81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 Set runtime options 83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 84 PetscCall(TSSetFromOptions(ts)); 85 86 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 87 Solve nonlinear system 88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 89 PetscCall(TSSolve(ts, u)); 90 PetscCall(TSGetSolveTime(ts, &ftime)); 91 PetscCall(TSGetStepNumber(ts, &steps)); 92 93 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 94 Free work space. 95 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 96 PetscCall(MatDestroy(&J)); 97 PetscCall(VecDestroy(&u)); 98 PetscCall(VecDestroy(&r)); 99 PetscCall(TSDestroy(&ts)); 100 PetscCall(DMDestroy(&da)); 101 102 PetscCall(PetscFinalize()); 103 return 0; 104 } 105 /* ------------------------------------------------------------------- */ 106 /* 107 RHSFunction - Evaluates nonlinear function, F(u). 108 109 Input Parameters: 110 . ts - the TS context 111 . U - input vector 112 . ptr - optional user-defined context, as set by TSSetFunction() 113 114 Output Parameter: 115 . F - function vector 116 */ 117 PetscErrorCode RHSFunction(TS ts, PetscReal ftime, Vec U, Vec F, void *ptr) 118 { 119 /* PETSC_UNUSED AppCtx *user=(AppCtx*)ptr; */ 120 DM da; 121 PetscInt i, j, Mx, My, xs, ys, xm, ym; 122 PetscReal two = 2.0, hx, hy, sx, sy; 123 PetscScalar u, uxx, uyy, **uarray, **f; 124 Vec localU; 125 126 PetscFunctionBeginUser; 127 PetscCall(TSGetDM(ts, &da)); 128 PetscCall(DMGetLocalVector(da, &localU)); 129 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)); 130 131 hx = 1.0 / (PetscReal)(Mx - 1); 132 sx = 1.0 / (hx * hx); 133 hy = 1.0 / (PetscReal)(My - 1); 134 sy = 1.0 / (hy * hy); 135 136 /* 137 Scatter ghost points to local vector,using the 2-step process 138 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 139 By placing code between these two statements, computations can be 140 done while messages are in transition. 141 */ 142 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU)); 143 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU)); 144 145 /* Get pointers to vector data */ 146 PetscCall(DMDAVecGetArrayRead(da, localU, &uarray)); 147 PetscCall(DMDAVecGetArray(da, F, &f)); 148 149 /* Get local grid boundaries */ 150 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 151 152 /* Compute function over the locally owned part of the grid */ 153 for (j = ys; j < ys + ym; j++) { 154 for (i = xs; i < xs + xm; i++) { 155 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) { 156 f[j][i] = uarray[j][i]; 157 continue; 158 } 159 u = uarray[j][i]; 160 uxx = (-two * u + uarray[j][i - 1] + uarray[j][i + 1]) * sx; 161 uyy = (-two * u + uarray[j - 1][i] + uarray[j + 1][i]) * sy; 162 f[j][i] = uxx + uyy; 163 } 164 } 165 166 /* Restore vectors */ 167 PetscCall(DMDAVecRestoreArrayRead(da, localU, &uarray)); 168 PetscCall(DMDAVecRestoreArray(da, F, &f)); 169 PetscCall(DMRestoreLocalVector(da, &localU)); 170 PetscCall(PetscLogFlops(11.0 * ym * xm)); 171 PetscFunctionReturn(PETSC_SUCCESS); 172 } 173 174 /* --------------------------------------------------------------------- */ 175 /* 176 RHSJacobian - User-provided routine to compute the Jacobian of 177 the nonlinear right-hand-side function of the ODE. 178 179 Input Parameters: 180 ts - the TS context 181 t - current time 182 U - global input vector 183 dummy - optional user-defined context, as set by TSetRHSJacobian() 184 185 Output Parameters: 186 J - Jacobian matrix 187 Jpre - optionally different preconditioning matrix 188 str - flag indicating matrix structure 189 */ 190 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx) 191 { 192 DM da; 193 DMDALocalInfo info; 194 PetscInt i, j; 195 PetscReal hx, hy, sx, sy; 196 197 PetscFunctionBeginUser; 198 PetscCall(TSGetDM(ts, &da)); 199 PetscCall(DMDAGetLocalInfo(da, &info)); 200 hx = 1.0 / (PetscReal)(info.mx - 1); 201 sx = 1.0 / (hx * hx); 202 hy = 1.0 / (PetscReal)(info.my - 1); 203 sy = 1.0 / (hy * hy); 204 for (j = info.ys; j < info.ys + info.ym; j++) { 205 for (i = info.xs; i < info.xs + info.xm; i++) { 206 PetscInt nc = 0; 207 MatStencil row, col[5]; 208 PetscScalar val[5]; 209 row.i = i; 210 row.j = j; 211 if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) { 212 col[nc].i = i; 213 col[nc].j = j; 214 val[nc++] = 1.0; 215 } else { 216 col[nc].i = i - 1; 217 col[nc].j = j; 218 val[nc++] = sx; 219 col[nc].i = i + 1; 220 col[nc].j = j; 221 val[nc++] = sx; 222 col[nc].i = i; 223 col[nc].j = j - 1; 224 val[nc++] = sy; 225 col[nc].i = i; 226 col[nc].j = j + 1; 227 val[nc++] = sy; 228 col[nc].i = i; 229 col[nc].j = j; 230 val[nc++] = -2 * sx - 2 * sy; 231 } 232 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES)); 233 } 234 } 235 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); 236 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); 237 if (J != Jpre) { 238 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); 239 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); 240 } 241 PetscFunctionReturn(PETSC_SUCCESS); 242 } 243 244 /* ------------------------------------------------------------------- */ 245 PetscErrorCode FormInitialSolution(DM da, Vec U, void *ptr) 246 { 247 AppCtx *user = (AppCtx *)ptr; 248 PetscReal c = user->c; 249 PetscInt i, j, xs, ys, xm, ym, Mx, My; 250 PetscScalar **u; 251 PetscReal hx, hy, x, y, r; 252 253 PetscFunctionBeginUser; 254 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)); 255 256 hx = 1.0 / (PetscReal)(Mx - 1); 257 hy = 1.0 / (PetscReal)(My - 1); 258 259 /* Get pointers to vector data */ 260 PetscCall(DMDAVecGetArray(da, U, &u)); 261 262 /* Get local grid boundaries */ 263 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 264 265 /* Compute function over the locally owned part of the grid */ 266 for (j = ys; j < ys + ym; j++) { 267 y = j * hy; 268 for (i = xs; i < xs + xm; i++) { 269 x = i * hx; 270 r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5)); 271 if (r < .125) u[j][i] = PetscExpReal(c * r * r * r); 272 else u[j][i] = 0.0; 273 } 274 } 275 276 /* Restore vectors */ 277 PetscCall(DMDAVecRestoreArray(da, U, &u)); 278 PetscFunctionReturn(PETSC_SUCCESS); 279 } 280 281 /*TEST 282 283 test: 284 args: -ts_max_steps 5 -ts_monitor 285 286 test: 287 suffix: 2 288 args: -ts_max_steps 5 -ts_monitor 289 290 test: 291 suffix: 3 292 args: -ts_max_steps 5 -snes_fd_color -ts_monitor 293 294 TEST*/ 295