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 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 DM da; 126 PetscInt i, j, Mx, My, xs, ys, xm, ym; 127 PetscReal two = 2.0, hx, hy, sx, sy; 128 PetscScalar u, uxx, uyy, **x, **f; 129 Vec localX; 130 131 PetscFunctionBeginUser; 132 PetscCall(TSGetDM(ts, &da)); 133 PetscCall(DMGetLocalVector(da, &localX)); 134 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)); 135 136 hx = 1.0 / (PetscReal)(Mx - 1); 137 sx = 1.0 / (hx * hx); 138 hy = 1.0 / (PetscReal)(My - 1); 139 sy = 1.0 / (hy * hy); 140 141 /* 142 Scatter ghost points to local vector,using the 2-step process 143 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 144 By placing code between these two statements, computations can be 145 done while messages are in transition. 146 */ 147 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 148 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 149 150 /* 151 Get pointers to vector data 152 */ 153 PetscCall(DMDAVecGetArrayRead(da, localX, &x)); 154 PetscCall(DMDAVecGetArray(da, F, &f)); 155 156 /* 157 Get local grid boundaries 158 */ 159 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 160 161 /* 162 Compute function over the locally owned part of the grid 163 */ 164 for (j = ys; j < ys + ym; j++) { 165 for (i = xs; i < xs + xm; i++) { 166 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) { 167 f[j][i] = x[j][i]; 168 continue; 169 } 170 u = x[j][i]; 171 uxx = (two * u - x[j][i - 1] - x[j][i + 1]) * sx; 172 uyy = (two * u - x[j - 1][i] - x[j + 1][i]) * sy; 173 /* f[j][i] = -(uxx + uyy); */ 174 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); 175 } 176 } 177 178 /* 179 Restore vectors 180 */ 181 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x)); 182 PetscCall(DMDAVecRestoreArray(da, F, &f)); 183 PetscCall(DMRestoreLocalVector(da, &localX)); 184 PetscCall(PetscLogFlops(11.0 * ym * xm)); 185 PetscFunctionReturn(0); 186 } 187 188 /* ------------------------------------------------------------------- */ 189 PetscErrorCode FormInitialSolution(DM da, Vec U) { 190 PetscInt i, j, xs, ys, xm, ym, Mx, My; 191 PetscScalar **u; 192 PetscReal hx, hy, x, y, r; 193 194 PetscFunctionBeginUser; 195 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)); 196 197 hx = 1.0 / (PetscReal)(Mx - 1); 198 hy = 1.0 / (PetscReal)(My - 1); 199 200 /* 201 Get pointers to vector data 202 */ 203 PetscCall(DMDAVecGetArray(da, U, &u)); 204 205 /* 206 Get local grid boundaries 207 */ 208 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 209 210 /* 211 Compute function over the locally owned part of the grid 212 */ 213 for (j = ys; j < ys + ym; j++) { 214 y = j * hy; 215 for (i = xs; i < xs + xm; i++) { 216 x = i * hx; 217 r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5)); 218 if (r < .125) u[j][i] = PetscExpReal(-30.0 * r * r * r); 219 else u[j][i] = 0.0; 220 } 221 } 222 223 /* 224 Restore vectors 225 */ 226 PetscCall(DMDAVecRestoreArray(da, U, &u)); 227 PetscFunctionReturn(0); 228 } 229 230 PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx) { 231 PetscReal norm; 232 MPI_Comm comm; 233 234 PetscFunctionBeginUser; 235 if (step < 0) PetscFunctionReturn(0); /* step of -1 indicates an interpolated solution */ 236 PetscCall(VecNorm(v, NORM_2, &norm)); 237 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm)); 238 PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm)); 239 PetscFunctionReturn(0); 240 } 241 242 /* 243 MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES. 244 Input Parameters: 245 snes - the SNES context 246 its - iteration number 247 fnorm - 2-norm function value (may be estimated) 248 ctx - optional user-defined context for private data for the 249 monitor routine, as set by SNESMonitorSet() 250 */ 251 PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf) { 252 PetscFunctionBeginUser; 253 PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf)); 254 PetscFunctionReturn(0); 255 } 256 257 /*TEST 258 259 test: 260 args: -ts_max_steps 5 261 262 test: 263 suffix: 2 264 args: -ts_max_steps 5 -snes_mf_operator 265 266 test: 267 suffix: 3 268 args: -ts_max_steps 5 -snes_mf -pc_type none 269 270 TEST*/ 271