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
main(int argc,char ** argv)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 */
RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void * ptr)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 matrix used to construct the preconditioner
186
187 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat J,Mat Jpre,PetscCtx ctx)188 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, PetscCtx 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
FormInitialSolution(DM da,Vec U,void * ptr)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