1 static char help[] = "Solves a time-dependent nonlinear PDE.\n";
2
3 /* ------------------------------------------------------------------------
4
5 This program solves the two-dimensional time-dependent Bratu problem
6 u_t = u_xx + u_yy + \lambda*exp(u),
7 on the domain 0 <= x,y <= 1,
8 with the boundary conditions
9 u(t,0,y) = 0, u_x(t,1,y) = 0,
10 u(t,x,0) = 0, u_x(t,x,1) = 0,
11 and the initial condition
12 u(0,x,y) = 0.
13 We discretize the right-hand side using finite differences with
14 uniform grid spacings hx,hy:
15 u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(hx^2)
16 u_yy = (u_{j+1} - 2u_{j} + u_{j-1})/(hy^2)
17
18 ------------------------------------------------------------------------- */
19
20 #include <petscdmda.h>
21 #include <petscts.h>
22
23 /*
24 User-defined application context - contains data needed by the
25 application-provided call-back routines.
26 */
27 typedef struct {
28 PetscReal lambda;
29 } AppCtx;
30
31 /*
32 FormIFunctionLocal - Evaluates nonlinear implicit function on local process patch
33 */
FormIFunctionLocal(DMDALocalInfo * info,PetscReal t,PetscScalar ** x,PetscScalar ** xdot,PetscScalar ** f,AppCtx * app)34 static PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal t, PetscScalar **x, PetscScalar **xdot, PetscScalar **f, AppCtx *app)
35 {
36 PetscInt i, j;
37 PetscReal lambda, hx, hy;
38 PetscScalar ut, u, ue, uw, un, us, uxx, uyy;
39
40 PetscFunctionBeginUser;
41 lambda = app->lambda;
42 hx = 1.0 / (PetscReal)(info->mx - 1);
43 hy = 1.0 / (PetscReal)(info->my - 1);
44
45 /*
46 Compute RHS function over the locally owned part of the grid
47 */
48 for (j = info->ys; j < info->ys + info->ym; j++) {
49 for (i = info->xs; i < info->xs + info->xm; i++) {
50 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
51 /* boundary points */
52 f[j][i] = x[j][i] - (PetscReal)0;
53 } else {
54 /* interior points */
55 ut = xdot[j][i];
56 u = x[j][i];
57 uw = x[j][i - 1];
58 ue = x[j][i + 1];
59 un = x[j + 1][i];
60 us = x[j - 1][i];
61
62 uxx = (uw - 2.0 * u + ue) / (hx * hx);
63 uyy = (un - 2.0 * u + us) / (hy * hy);
64 f[j][i] = ut - uxx - uyy - lambda * PetscExpScalar(u);
65 }
66 }
67 }
68 PetscFunctionReturn(PETSC_SUCCESS);
69 }
70
71 /*
72 FormIJacobianLocal - Evaluates implicit Jacobian matrix on local process patch
73 */
FormIJacobianLocal(DMDALocalInfo * info,PetscReal t,PetscScalar ** x,PetscScalar ** xdot,PetscScalar shift,Mat jac,Mat jacpre,AppCtx * app)74 static PetscErrorCode FormIJacobianLocal(DMDALocalInfo *info, PetscReal t, PetscScalar **x, PetscScalar **xdot, PetscScalar shift, Mat jac, Mat jacpre, AppCtx *app)
75 {
76 PetscInt i, j, k;
77 MatStencil col[5], row;
78 PetscScalar v[5], lambda, hx, hy;
79
80 PetscFunctionBeginUser;
81 lambda = app->lambda;
82 hx = 1.0 / (PetscReal)(info->mx - 1);
83 hy = 1.0 / (PetscReal)(info->my - 1);
84
85 /*
86 Compute Jacobian entries for the locally owned part of the grid
87 */
88 for (j = info->ys; j < info->ys + info->ym; j++) {
89 for (i = info->xs; i < info->xs + info->xm; i++) {
90 row.j = j;
91 row.i = i;
92 k = 0;
93 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) {
94 /* boundary points */
95 v[0] = 1.0;
96 PetscCall(MatSetValuesStencil(jacpre, 1, &row, 1, &row, v, INSERT_VALUES));
97 } else {
98 /* interior points */
99 v[k] = -1.0 / (hy * hy);
100 col[k].j = j - 1;
101 col[k].i = i;
102 k++;
103 v[k] = -1.0 / (hx * hx);
104 col[k].j = j;
105 col[k].i = i - 1;
106 k++;
107
108 v[k] = shift + 2.0 / (hx * hx) + 2.0 / (hy * hy) - lambda * PetscExpScalar(x[j][i]);
109 col[k].j = j;
110 col[k].i = i;
111 k++;
112
113 v[k] = -1.0 / (hx * hx);
114 col[k].j = j;
115 col[k].i = i + 1;
116 k++;
117 v[k] = -1.0 / (hy * hy);
118 col[k].j = j + 1;
119 col[k].i = i;
120 k++;
121
122 PetscCall(MatSetValuesStencil(jacpre, 1, &row, k, col, v, INSERT_VALUES));
123 }
124 }
125 }
126
127 /*
128 Assemble matrix
129 */
130 PetscCall(MatAssemblyBegin(jacpre, MAT_FINAL_ASSEMBLY));
131 PetscCall(MatAssemblyEnd(jacpre, MAT_FINAL_ASSEMBLY));
132 PetscFunctionReturn(PETSC_SUCCESS);
133 }
134
main(int argc,char ** argv)135 int main(int argc, char **argv)
136 {
137 TS ts; /* ODE integrator */
138 DM da; /* DM context */
139 Vec U; /* solution vector */
140 DMBoundaryType bt = DM_BOUNDARY_NONE;
141 DMDAStencilType st = DMDA_STENCIL_STAR;
142 PetscInt sw = 1;
143 PetscInt N = 17;
144 PetscInt n = PETSC_DECIDE;
145 AppCtx app;
146
147 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148 Initialize program
149 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150 PetscFunctionBeginUser;
151 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
152 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex21 options", "");
153 {
154 app.lambda = 6.8;
155 app.lambda = 6.0;
156 PetscCall(PetscOptionsReal("-lambda", "", "", app.lambda, &app.lambda, NULL));
157 }
158 PetscOptionsEnd();
159
160 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161 Create DM context
162 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bt, bt, st, N, N, n, n, 1, sw, NULL, NULL, &da));
164 PetscCall(DMSetFromOptions(da));
165 PetscCall(DMSetUp(da));
166 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0, 1.0));
167
168 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169 Create timestepping solver context
170 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
172 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
173 PetscCall(TSSetDM(ts, da));
174 PetscCall(DMDestroy(&da));
175
176 PetscCall(TSGetDM(ts, &da));
177 PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (DMDATSIFunctionLocalFn *)FormIFunctionLocal, &app));
178 PetscCall(DMDATSSetIJacobianLocal(da, (DMDATSIJacobianLocalFn *)FormIJacobianLocal, &app));
179
180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181 Set solver options
182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183 PetscCall(TSSetType(ts, TSBDF));
184 PetscCall(TSSetTimeStep(ts, 1e-4));
185 PetscCall(TSSetMaxTime(ts, 1.0));
186 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
187 PetscCall(TSSetFromOptions(ts));
188
189 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190 Set initial conditions
191 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192 PetscCall(TSGetDM(ts, &da));
193 PetscCall(DMCreateGlobalVector(da, &U));
194 PetscCall(VecSet(U, 0.0));
195 PetscCall(TSSetSolution(ts, U));
196
197 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198 Run timestepping solver
199 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200 PetscCall(TSSolve(ts, U));
201
202 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203 All PETSc objects should be destroyed when they are no longer needed.
204 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205 PetscCall(VecDestroy(&U));
206 PetscCall(TSDestroy(&ts));
207 PetscCall(PetscFinalize());
208 return 0;
209 }
210
211 /*TEST
212
213 testset:
214 requires: !single !complex
215 args: -da_grid_x 5 -da_grid_y 5 -da_refine 2 -dm_view -ts_type bdf -ts_adapt_type none -ts_time_step 1e-3 -ts_monitor -ts_max_steps 5 -ts_view -snes_rtol 1e-6 -snes_type ngmres -npc_snes_type fas
216 filter: grep -v "total number of"
217 test:
218 suffix: 1_bdf_ngmres_fas_ms
219 args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop
220 test:
221 suffix: 2_bdf_ngmres_fas_ms
222 args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop
223 nsize: 2
224 test:
225 suffix: 1_bdf_ngmres_fas_ngs
226 args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop
227 test:
228 suffix: 2_bdf_ngmres_fas_ngs
229 args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop
230 nsize: 2
231
232 TEST*/
233