Actual source code: ex21.c
1: static char help[] = "Solves a time-dependent nonlinear PDE.\n";
3: /* ------------------------------------------------------------------------
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)
18: ------------------------------------------------------------------------- */
20: #include <petscdmda.h>
21: #include <petscts.h>
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;
31: /*
32: FormIFunctionLocal - Evaluates nonlinear implicit function on local process patch
33: */
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;
40: PetscFunctionBeginUser;
41: lambda = app->lambda;
42: hx = 1.0 / (PetscReal)(info->mx - 1);
43: hy = 1.0 / (PetscReal)(info->my - 1);
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];
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: }
71: /*
72: FormIJacobianLocal - Evaluates implicit Jacobian matrix on local process patch
73: */
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;
80: PetscFunctionBeginUser;
81: lambda = app->lambda;
82: hx = 1.0 / (PetscReal)(info->mx - 1);
83: hy = 1.0 / (PetscReal)(info->my - 1);
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++;
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++;
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++;
122: PetscCall(MatSetValuesStencil(jacpre, 1, &row, k, col, v, INSERT_VALUES));
123: }
124: }
125: }
127: /*
128: Assemble matrix
129: */
130: PetscCall(MatAssemblyBegin(jacpre, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatAssemblyEnd(jacpre, MAT_FINAL_ASSEMBLY));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
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;
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();
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));
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));
176: PetscCall(TSGetDM(ts, &da));
177: PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (DMDATSIFunctionLocalFn *)FormIFunctionLocal, &app));
178: PetscCall(DMDATSSetIJacobianLocal(da, (DMDATSIJacobianLocalFn *)FormIJacobianLocal, &app));
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));
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));
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Run timestepping solver
199: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200: PetscCall(TSSolve(ts, U));
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: }
211: /*TEST
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_dt 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
232: TEST*/