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 PetscInt i, j; 36 PetscReal lambda, hx, hy; 37 PetscScalar ut, u, ue, uw, un, us, uxx, uyy; 38 39 PetscFunctionBeginUser; 40 lambda = app->lambda; 41 hx = 1.0 / (PetscReal)(info->mx - 1); 42 hy = 1.0 / (PetscReal)(info->my - 1); 43 44 /* 45 Compute RHS function over the locally owned part of the grid 46 */ 47 for (j = info->ys; j < info->ys + info->ym; j++) { 48 for (i = info->xs; i < info->xs + info->xm; i++) { 49 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) { 50 /* boundary points */ 51 f[j][i] = x[j][i] - (PetscReal)0; 52 } else { 53 /* interior points */ 54 ut = xdot[j][i]; 55 u = x[j][i]; 56 uw = x[j][i - 1]; 57 ue = x[j][i + 1]; 58 un = x[j + 1][i]; 59 us = x[j - 1][i]; 60 61 uxx = (uw - 2.0 * u + ue) / (hx * hx); 62 uyy = (un - 2.0 * u + us) / (hy * hy); 63 f[j][i] = ut - uxx - uyy - lambda * PetscExpScalar(u); 64 } 65 } 66 } 67 PetscFunctionReturn(0); 68 } 69 70 /* 71 FormIJacobianLocal - Evaluates implicit Jacobian matrix on local process patch 72 */ 73 static PetscErrorCode FormIJacobianLocal(DMDALocalInfo *info, PetscReal t, PetscScalar **x, PetscScalar **xdot, PetscScalar shift, Mat jac, Mat jacpre, AppCtx *app) { 74 PetscInt i, j, k; 75 MatStencil col[5], row; 76 PetscScalar v[5], lambda, hx, hy; 77 78 PetscFunctionBeginUser; 79 lambda = app->lambda; 80 hx = 1.0 / (PetscReal)(info->mx - 1); 81 hy = 1.0 / (PetscReal)(info->my - 1); 82 83 /* 84 Compute Jacobian entries for the locally owned part of the grid 85 */ 86 for (j = info->ys; j < info->ys + info->ym; j++) { 87 for (i = info->xs; i < info->xs + info->xm; i++) { 88 row.j = j; 89 row.i = i; 90 k = 0; 91 if (i == 0 || j == 0 || i == info->mx - 1 || j == info->my - 1) { 92 /* boundary points */ 93 v[0] = 1.0; 94 PetscCall(MatSetValuesStencil(jacpre, 1, &row, 1, &row, v, INSERT_VALUES)); 95 } else { 96 /* interior points */ 97 v[k] = -1.0 / (hy * hy); 98 col[k].j = j - 1; 99 col[k].i = i; 100 k++; 101 v[k] = -1.0 / (hx * hx); 102 col[k].j = j; 103 col[k].i = i - 1; 104 k++; 105 106 v[k] = shift + 2.0 / (hx * hx) + 2.0 / (hy * hy) - lambda * PetscExpScalar(x[j][i]); 107 col[k].j = j; 108 col[k].i = i; 109 k++; 110 111 v[k] = -1.0 / (hx * hx); 112 col[k].j = j; 113 col[k].i = i + 1; 114 k++; 115 v[k] = -1.0 / (hy * hy); 116 col[k].j = j + 1; 117 col[k].i = i; 118 k++; 119 120 PetscCall(MatSetValuesStencil(jacpre, 1, &row, k, col, v, INSERT_VALUES)); 121 } 122 } 123 } 124 125 /* 126 Assemble matrix 127 */ 128 PetscCall(MatAssemblyBegin(jacpre, MAT_FINAL_ASSEMBLY)); 129 PetscCall(MatAssemblyEnd(jacpre, MAT_FINAL_ASSEMBLY)); 130 PetscFunctionReturn(0); 131 } 132 133 int main(int argc, char **argv) { 134 TS ts; /* ODE integrator */ 135 DM da; /* DM context */ 136 Vec U; /* solution vector */ 137 DMBoundaryType bt = DM_BOUNDARY_NONE; 138 DMDAStencilType st = DMDA_STENCIL_STAR; 139 PetscInt sw = 1; 140 PetscInt N = 17; 141 PetscInt n = PETSC_DECIDE; 142 AppCtx app; 143 144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 Initialize program 146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 147 PetscFunctionBeginUser; 148 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 149 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex21 options", ""); 150 { 151 app.lambda = 6.8; 152 app.lambda = 6.0; 153 PetscCall(PetscOptionsReal("-lambda", "", "", app.lambda, &app.lambda, NULL)); 154 } 155 PetscOptionsEnd(); 156 157 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 158 Create DM context 159 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 160 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, bt, bt, st, N, N, n, n, 1, sw, NULL, NULL, &da)); 161 PetscCall(DMSetFromOptions(da)); 162 PetscCall(DMSetUp(da)); 163 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0, 1.0)); 164 165 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 166 Create timestepping solver context 167 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 168 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 169 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 170 PetscCall(TSSetDM(ts, da)); 171 PetscCall(DMDestroy(&da)); 172 173 PetscCall(TSGetDM(ts, &da)); 174 PetscCall(DMDATSSetIFunctionLocal(da, INSERT_VALUES, (DMDATSIFunctionLocal)FormIFunctionLocal, &app)); 175 PetscCall(DMDATSSetIJacobianLocal(da, (DMDATSIJacobianLocal)FormIJacobianLocal, &app)); 176 177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178 Set solver options 179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 180 PetscCall(TSSetType(ts, TSBDF)); 181 PetscCall(TSSetTimeStep(ts, 1e-4)); 182 PetscCall(TSSetMaxTime(ts, 1.0)); 183 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 184 PetscCall(TSSetFromOptions(ts)); 185 186 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 187 Set initial conditions 188 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 189 PetscCall(TSGetDM(ts, &da)); 190 PetscCall(DMCreateGlobalVector(da, &U)); 191 PetscCall(VecSet(U, 0.0)); 192 PetscCall(TSSetSolution(ts, U)); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Run timestepping solver 196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197 PetscCall(TSSolve(ts, U)); 198 199 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 200 All PETSc objects should be destroyed when they are no longer needed. 201 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 202 PetscCall(VecDestroy(&U)); 203 PetscCall(TSDestroy(&ts)); 204 PetscCall(PetscFinalize()); 205 return 0; 206 } 207 208 /*TEST 209 210 testset: 211 requires: !single !complex 212 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 213 filter: grep -v "total number of" 214 test: 215 suffix: 1_bdf_ngmres_fas_ms 216 args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop 217 test: 218 suffix: 2_bdf_ngmres_fas_ms 219 args: -prefix_push npc_fas_levels_ -snes_type ms -snes_max_it 5 -ksp_type preonly -prefix_pop 220 nsize: 2 221 test: 222 suffix: 1_bdf_ngmres_fas_ngs 223 args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop 224 test: 225 suffix: 2_bdf_ngmres_fas_ngs 226 args: -prefix_push npc_fas_levels_ -snes_type ngs -snes_max_it 5 -prefix_pop 227 nsize: 2 228 229 TEST*/ 230