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