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