1 static char help[] = "Nonlinear, time-dependent PDE in 2d.\n"; 2 3 /* 4 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 5 Include "petscts.h" so that we can use SNES solvers. Note that this 6 file automatically includes: 7 petscsys.h - base PETSc routines petscvec.h - vectors 8 petscmat.h - matrices 9 petscis.h - index sets petscksp.h - Krylov subspace methods 10 petscviewer.h - viewers petscpc.h - preconditioners 11 petscksp.h - linear solvers 12 */ 13 #include <petscdm.h> 14 #include <petscdmda.h> 15 #include <petscts.h> 16 17 /* 18 User-defined routines 19 */ 20 extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec); 21 extern PetscErrorCode MyTSMonitor(TS, PetscInt, PetscReal, Vec, void *); 22 extern PetscErrorCode MySNESMonitor(SNES, PetscInt, PetscReal, PetscViewerAndFormat *); 23 24 int main(int argc, char **argv) 25 { 26 TS ts; /* time integrator */ 27 SNES snes; 28 Vec x, r; /* solution, residual vectors */ 29 DM da; 30 PetscViewerAndFormat *vf; 31 32 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 33 Initialize program 34 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 35 PetscFunctionBeginUser; 36 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 37 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 38 Create distributed array (DMDA) to manage parallel grid and vectors 39 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 40 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)); 41 PetscCall(DMSetFromOptions(da)); 42 PetscCall(DMSetUp(da)); 43 44 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 45 Extract global vectors from DMDA; then duplicate for remaining 46 vectors that are the same types 47 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 48 PetscCall(DMCreateGlobalVector(da, &x)); 49 PetscCall(VecDuplicate(x, &r)); 50 51 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 52 Create timestepping solver context 53 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 54 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 55 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 56 PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, da)); 57 58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 59 Create matrix data structure; set Jacobian evaluation routine 60 61 Set Jacobian matrix data structure and default Jacobian evaluation 62 routine. User can override with: 63 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 64 (unless user explicitly sets preconditioner) 65 -snes_mf_operator : form preconditioning matrix as set by the user, 66 but use matrix-free approx for Jacobian-vector 67 products within Newton-Krylov method 68 69 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 70 71 PetscCall(TSSetMaxTime(ts, 1.0)); 72 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 73 PetscCall(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL)); 74 PetscCall(TSSetDM(ts, da)); 75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 76 Customize nonlinear solver 77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 78 PetscCall(TSSetType(ts, TSBEULER)); 79 PetscCall(TSGetSNES(ts, &snes)); 80 PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf)); 81 PetscCall(SNESMonitorSet(snes, (PetscErrorCode(*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy)); 82 83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 84 Set initial conditions 85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 86 PetscCall(FormInitialSolution(da, x)); 87 PetscCall(TSSetTimeStep(ts, .0001)); 88 PetscCall(TSSetSolution(ts, x)); 89 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Set runtime options 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 93 PetscCall(TSSetFromOptions(ts)); 94 95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 96 Solve nonlinear system 97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 98 PetscCall(TSSolve(ts, x)); 99 100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 Free work space. All PETSc objects should be destroyed when they 102 are no longer needed. 103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 104 PetscCall(VecDestroy(&x)); 105 PetscCall(VecDestroy(&r)); 106 PetscCall(TSDestroy(&ts)); 107 PetscCall(DMDestroy(&da)); 108 109 PetscCall(PetscFinalize()); 110 return 0; 111 } 112 /* ------------------------------------------------------------------- */ 113 /* 114 FormFunction - Evaluates nonlinear function, F(x). 115 116 Input Parameters: 117 . ts - the TS context 118 . X - input vector 119 . ptr - optional user-defined context, as set by SNESSetFunction() 120 121 Output Parameter: 122 . F - function vector 123 */ 124 PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr) 125 { 126 DM da; 127 PetscInt i, j, Mx, My, xs, ys, xm, ym; 128 PetscReal two = 2.0, hx, hy, sx, sy; 129 PetscScalar u, uxx, uyy, **x, **f; 130 Vec localX; 131 132 PetscFunctionBeginUser; 133 PetscCall(TSGetDM(ts, &da)); 134 PetscCall(DMGetLocalVector(da, &localX)); 135 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)); 136 137 hx = 1.0 / (PetscReal)(Mx - 1); 138 sx = 1.0 / (hx * hx); 139 hy = 1.0 / (PetscReal)(My - 1); 140 sy = 1.0 / (hy * hy); 141 142 /* 143 Scatter ghost points to local vector,using the 2-step process 144 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 145 By placing code between these two statements, computations can be 146 done while messages are in transition. 147 */ 148 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 149 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 150 151 /* 152 Get pointers to vector data 153 */ 154 PetscCall(DMDAVecGetArrayRead(da, localX, &x)); 155 PetscCall(DMDAVecGetArray(da, F, &f)); 156 157 /* 158 Get local grid boundaries 159 */ 160 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 161 162 /* 163 Compute function over the locally owned part of the grid 164 */ 165 for (j = ys; j < ys + ym; j++) { 166 for (i = xs; i < xs + xm; i++) { 167 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) { 168 f[j][i] = x[j][i]; 169 continue; 170 } 171 u = x[j][i]; 172 uxx = (two * u - x[j][i - 1] - x[j][i + 1]) * sx; 173 uyy = (two * u - x[j - 1][i] - x[j + 1][i]) * sy; 174 /* f[j][i] = -(uxx + uyy); */ 175 f[j][i] = -u * (uxx + uyy) - (4.0 - 1.0) * ((x[j][i + 1] - x[j][i - 1]) * (x[j][i + 1] - x[j][i - 1]) * .25 * sx + (x[j + 1][i] - x[j - 1][i]) * (x[j + 1][i] - x[j - 1][i]) * .25 * sy); 176 } 177 } 178 179 /* 180 Restore vectors 181 */ 182 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x)); 183 PetscCall(DMDAVecRestoreArray(da, F, &f)); 184 PetscCall(DMRestoreLocalVector(da, &localX)); 185 PetscCall(PetscLogFlops(11.0 * ym * xm)); 186 PetscFunctionReturn(PETSC_SUCCESS); 187 } 188 189 /* ------------------------------------------------------------------- */ 190 PetscErrorCode FormInitialSolution(DM da, Vec U) 191 { 192 PetscInt i, j, xs, ys, xm, ym, Mx, My; 193 PetscScalar **u; 194 PetscReal hx, hy, x, y, r; 195 196 PetscFunctionBeginUser; 197 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)); 198 199 hx = 1.0 / (PetscReal)(Mx - 1); 200 hy = 1.0 / (PetscReal)(My - 1); 201 202 /* 203 Get pointers to vector data 204 */ 205 PetscCall(DMDAVecGetArray(da, U, &u)); 206 207 /* 208 Get local grid boundaries 209 */ 210 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 211 212 /* 213 Compute function over the locally owned part of the grid 214 */ 215 for (j = ys; j < ys + ym; j++) { 216 y = j * hy; 217 for (i = xs; i < xs + xm; i++) { 218 x = i * hx; 219 r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5)); 220 if (r < .125) u[j][i] = PetscExpReal(-30.0 * r * r * r); 221 else u[j][i] = 0.0; 222 } 223 } 224 225 /* 226 Restore vectors 227 */ 228 PetscCall(DMDAVecRestoreArray(da, U, &u)); 229 PetscFunctionReturn(PETSC_SUCCESS); 230 } 231 232 PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx) 233 { 234 PetscReal norm; 235 MPI_Comm comm; 236 237 PetscFunctionBeginUser; 238 if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* step of -1 indicates an interpolated solution */ 239 PetscCall(VecNorm(v, NORM_2, &norm)); 240 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm)); 241 PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm)); 242 PetscFunctionReturn(PETSC_SUCCESS); 243 } 244 245 /* 246 MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES. 247 Input Parameters: 248 snes - the SNES context 249 its - iteration number 250 fnorm - 2-norm function value (may be estimated) 251 ctx - optional user-defined context for private data for the 252 monitor routine, as set by SNESMonitorSet() 253 */ 254 PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf) 255 { 256 PetscFunctionBeginUser; 257 PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf)); 258 PetscFunctionReturn(PETSC_SUCCESS); 259 } 260 261 /*TEST 262 263 test: 264 args: -ts_max_steps 5 265 266 test: 267 suffix: 2 268 args: -ts_max_steps 5 -snes_mf_operator 269 270 test: 271 suffix: 3 272 args: -ts_max_steps 5 -snes_mf -pc_type none 273 274 TEST*/ 275