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