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