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 */ 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 */ 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 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