1 static char help[] = "Solves heat equation in 1d.\n"; 2 3 /* 4 Solves the equation 5 6 u_t = kappa \Delta u 7 Periodic boundary conditions 8 9 Evolve the heat equation: 10 --------------- 11 ./heat -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor 12 13 Evolve the Allen-Cahn equation: 14 --------------- 15 ./heat -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -ts_type cn -da_refine 5 -allen-cahn -kappa .001 -ts_max_time 5 -mymonitor 16 17 Evolve the Allen-Cahn equation: zoom in on part of the domain 18 --------------- 19 ./heat -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -allen-cahn -kappa .001 -ts_max_time 5 -zoom .25,.45 -mymonitor 20 21 The option -square_initial indicates it should use a square wave initial condition otherwise it loads the file InitialSolution.heat as the initial solution. You should run with 22 ./heat -square_initial -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 9 -ts_max_time 1.e-4 -ts_dt .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25 -ts_max_steps 15 23 to generate InitialSolution.heat 24 25 */ 26 #include <petscdm.h> 27 #include <petscdmda.h> 28 #include <petscts.h> 29 #include <petscdraw.h> 30 31 /* 32 User-defined routines 33 */ 34 extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, void *), MyDestroy(void **); 35 typedef struct { 36 PetscReal kappa; 37 PetscBool allencahn; 38 PetscDrawViewPorts *ports; 39 } UserCtx; 40 41 int main(int argc, char **argv) 42 { 43 TS ts; /* time integrator */ 44 Vec x, r; /* solution, residual vectors */ 45 PetscInt steps, Mx; 46 DM da; 47 PetscReal dt; 48 UserCtx ctx; 49 PetscBool mymonitor; 50 PetscViewer viewer; 51 PetscBool flg; 52 53 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 54 Initialize program 55 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 56 PetscFunctionBeginUser; 57 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 58 ctx.kappa = 1.0; 59 PetscCall(PetscOptionsGetReal(NULL, NULL, "-kappa", &ctx.kappa, NULL)); 60 ctx.allencahn = PETSC_FALSE; 61 PetscCall(PetscOptionsHasName(NULL, NULL, "-allen-cahn", &ctx.allencahn)); 62 PetscCall(PetscOptionsHasName(NULL, NULL, "-mymonitor", &mymonitor)); 63 64 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 65 Create distributed array (DMDA) to manage parallel grid and vectors 66 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 67 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10, 1, 2, NULL, &da)); 68 PetscCall(DMSetFromOptions(da)); 69 PetscCall(DMSetUp(da)); 70 PetscCall(DMDASetFieldName(da, 0, "Heat equation: u")); 71 PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 72 dt = 1.0 / (ctx.kappa * Mx * Mx); 73 74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 75 Extract global vectors from DMDA; then duplicate for remaining 76 vectors that are the same types 77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 78 PetscCall(DMCreateGlobalVector(da, &x)); 79 PetscCall(VecDuplicate(x, &r)); 80 81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 Create timestepping solver context 83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 84 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 85 PetscCall(TSSetDM(ts, da)); 86 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 87 PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, &ctx)); 88 89 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 90 Customize nonlinear solver 91 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 92 PetscCall(TSSetType(ts, TSCN)); 93 94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 95 Set initial conditions 96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 97 PetscCall(FormInitialSolution(da, x)); 98 PetscCall(TSSetTimeStep(ts, dt)); 99 PetscCall(TSSetMaxTime(ts, .02)); 100 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE)); 101 PetscCall(TSSetSolution(ts, x)); 102 103 if (mymonitor) { 104 ctx.ports = NULL; 105 PetscCall(TSMonitorSet(ts, MyMonitor, &ctx, MyDestroy)); 106 } 107 108 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 109 Set runtime options 110 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 111 PetscCall(TSSetFromOptions(ts)); 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Solve nonlinear system 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 PetscCall(TSSolve(ts, x)); 117 PetscCall(TSGetStepNumber(ts, &steps)); 118 PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg)); 119 if (flg) { 120 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_WRITE, &viewer)); 121 PetscCall(VecView(x, viewer)); 122 PetscCall(PetscViewerDestroy(&viewer)); 123 } 124 125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 126 Free work space. All PETSc objects should be destroyed when they 127 are no longer needed. 128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 129 PetscCall(VecDestroy(&x)); 130 PetscCall(VecDestroy(&r)); 131 PetscCall(TSDestroy(&ts)); 132 PetscCall(DMDestroy(&da)); 133 134 PetscCall(PetscFinalize()); 135 return 0; 136 } 137 /* ------------------------------------------------------------------- */ 138 /* 139 FormFunction - Evaluates nonlinear function, F(x). 140 141 Input Parameters: 142 . ts - the TS context 143 . X - input vector 144 . ptr - optional user-defined context, as set by SNESSetFunction() 145 146 Output Parameter: 147 . F - function vector 148 */ 149 PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr) 150 { 151 DM da; 152 PetscInt i, Mx, xs, xm; 153 PetscReal hx, sx; 154 PetscScalar *x, *f; 155 Vec localX; 156 UserCtx *ctx = (UserCtx *)ptr; 157 158 PetscFunctionBegin; 159 PetscCall(TSGetDM(ts, &da)); 160 PetscCall(DMGetLocalVector(da, &localX)); 161 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 162 163 hx = 1.0 / (PetscReal)Mx; 164 sx = 1.0 / (hx * hx); 165 166 /* 167 Scatter ghost points to local vector,using the 2-step process 168 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 169 By placing code between these two statements, computations can be 170 done while messages are in transition. 171 */ 172 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 173 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 174 175 /* 176 Get pointers to vector data 177 */ 178 PetscCall(DMDAVecGetArrayRead(da, localX, &x)); 179 PetscCall(DMDAVecGetArray(da, F, &f)); 180 181 /* 182 Get local grid boundaries 183 */ 184 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 185 186 /* 187 Compute function over the locally owned part of the grid 188 */ 189 for (i = xs; i < xs + xm; i++) { 190 f[i] = ctx->kappa * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx; 191 if (ctx->allencahn) f[i] += (x[i] - x[i] * x[i] * x[i]); 192 } 193 194 /* 195 Restore vectors 196 */ 197 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x)); 198 PetscCall(DMDAVecRestoreArray(da, F, &f)); 199 PetscCall(DMRestoreLocalVector(da, &localX)); 200 PetscFunctionReturn(PETSC_SUCCESS); 201 } 202 203 /* ------------------------------------------------------------------- */ 204 PetscErrorCode FormInitialSolution(DM da, Vec U) 205 { 206 PetscInt i, xs, xm, Mx, scale = 1, N; 207 PetscScalar *u; 208 const PetscScalar *f; 209 PetscReal hx, x, r; 210 Vec finesolution; 211 PetscViewer viewer; 212 PetscBool flg; 213 214 PetscFunctionBegin; 215 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 216 217 hx = 1.0 / (PetscReal)Mx; 218 219 /* 220 Get pointers to vector data 221 */ 222 PetscCall(DMDAVecGetArray(da, U, &u)); 223 224 /* 225 Get local grid boundaries 226 */ 227 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 228 229 /* InitialSolution is obtained with 230 ./heat -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 9 -ts_max_time 1.e-4 -ts_dt .125e-6 -snes_atol 1.e-25 -snes_rtol 1.e-25 -ts_max_steps 15 231 */ 232 PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg)); 233 if (!flg) { 234 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer)); 235 PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution)); 236 PetscCall(VecLoad(finesolution, viewer)); 237 PetscCall(PetscViewerDestroy(&viewer)); 238 PetscCall(VecGetSize(finesolution, &N)); 239 scale = N / Mx; 240 PetscCall(VecGetArrayRead(finesolution, &f)); 241 } 242 243 /* 244 Compute function over the locally owned part of the grid 245 */ 246 for (i = xs; i < xs + xm; i++) { 247 x = i * hx; 248 r = PetscSqrtReal((x - .5) * (x - .5)); 249 if (r < .125) u[i] = 1.0; 250 else u[i] = -.5; 251 252 /* With the initial condition above the method is first order in space */ 253 /* this is a smooth initial condition so the method becomes second order in space */ 254 /*u[i] = PetscSinScalar(2*PETSC_PI*x); */ 255 /* u[i] = f[scale*i];*/ 256 if (!flg) u[i] = f[scale * i]; 257 } 258 if (!flg) { 259 PetscCall(VecRestoreArrayRead(finesolution, &f)); 260 PetscCall(VecDestroy(&finesolution)); 261 } 262 263 /* 264 Restore vectors 265 */ 266 PetscCall(DMDAVecRestoreArray(da, U, &u)); 267 PetscFunctionReturn(PETSC_SUCCESS); 268 } 269 270 /* 271 This routine is not parallel 272 */ 273 PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr) 274 { 275 UserCtx *ctx = (UserCtx *)ptr; 276 PetscDrawLG lg; 277 PetscScalar *u; 278 PetscInt Mx, i, xs, xm, cnt; 279 PetscReal x, y, hx, pause, sx, len, max, xx[2], yy[2]; 280 PetscDraw draw; 281 Vec localU; 282 DM da; 283 int colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE}; 284 const char *const legend[] = {"-kappa (\\grad u,\\grad u)", "(1 - u^2)^2"}; 285 PetscDrawAxis axis; 286 PetscDrawViewPorts *ports; 287 PetscReal vbounds[] = {-1.1, 1.1}; 288 289 PetscFunctionBegin; 290 PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds)); 291 PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1200, 800)); 292 PetscCall(TSGetDM(ts, &da)); 293 PetscCall(DMGetLocalVector(da, &localU)); 294 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 295 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 296 hx = 1.0 / (PetscReal)Mx; 297 sx = 1.0 / (hx * hx); 298 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU)); 299 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU)); 300 PetscCall(DMDAVecGetArrayRead(da, localU, &u)); 301 302 PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg)); 303 PetscCall(PetscDrawLGGetDraw(lg, &draw)); 304 PetscCall(PetscDrawCheckResizedWindow(draw)); 305 if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports)); 306 ports = ctx->ports; 307 PetscCall(PetscDrawLGGetAxis(lg, &axis)); 308 PetscCall(PetscDrawLGReset(lg)); 309 310 xx[0] = 0.0; 311 xx[1] = 1.0; 312 cnt = 2; 313 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL)); 314 xs = xx[0] / hx; 315 xm = (xx[1] - xx[0]) / hx; 316 317 /* 318 Plot the energies 319 */ 320 PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->allencahn ? 1 : 0))); 321 PetscCall(PetscDrawLGSetColors(lg, colors + 1)); 322 PetscCall(PetscDrawViewPortsSet(ports, 2)); 323 x = hx * xs; 324 for (i = xs; i < xs + xm; i++) { 325 xx[0] = xx[1] = x; 326 yy[0] = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx); 327 if (ctx->allencahn) yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i])); 328 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 329 x += hx; 330 } 331 PetscCall(PetscDrawGetPause(draw, &pause)); 332 PetscCall(PetscDrawSetPause(draw, 0.0)); 333 PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", "")); 334 PetscCall(PetscDrawLGSetLegend(lg, legend)); 335 PetscCall(PetscDrawLGDraw(lg)); 336 337 /* 338 Plot the forces 339 */ 340 PetscCall(PetscDrawViewPortsSet(ports, 1)); 341 PetscCall(PetscDrawLGReset(lg)); 342 x = xs * hx; 343 max = 0.; 344 for (i = xs; i < xs + xm; i++) { 345 xx[0] = xx[1] = x; 346 yy[0] = PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx); 347 max = PetscMax(max, PetscAbs(yy[0])); 348 if (ctx->allencahn) { 349 yy[1] = PetscRealPart(u[i] - u[i] * u[i] * u[i]); 350 max = PetscMax(max, PetscAbs(yy[1])); 351 } 352 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 353 x += hx; 354 } 355 PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", "")); 356 PetscCall(PetscDrawLGSetLegend(lg, NULL)); 357 PetscCall(PetscDrawLGDraw(lg)); 358 359 /* 360 Plot the solution 361 */ 362 PetscCall(PetscDrawLGSetDimension(lg, 1)); 363 PetscCall(PetscDrawViewPortsSet(ports, 0)); 364 PetscCall(PetscDrawLGReset(lg)); 365 x = hx * xs; 366 PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1)); 367 PetscCall(PetscDrawLGSetColors(lg, colors)); 368 for (i = xs; i < xs + xm; i++) { 369 xx[0] = x; 370 yy[0] = PetscRealPart(u[i]); 371 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 372 x += hx; 373 } 374 PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", "")); 375 PetscCall(PetscDrawLGDraw(lg)); 376 377 /* 378 Print the forces as arrows on the solution 379 */ 380 x = hx * xs; 381 cnt = xm / 60; 382 cnt = (!cnt) ? 1 : cnt; 383 384 for (i = xs; i < xs + xm; i += cnt) { 385 y = PetscRealPart(u[i]); 386 len = .5 * PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max; 387 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED)); 388 if (ctx->allencahn) { 389 len = .5 * PetscRealPart(u[i] - u[i] * u[i] * u[i]) / max; 390 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE)); 391 } 392 x += cnt * hx; 393 } 394 PetscCall(DMDAVecRestoreArrayRead(da, localU, &x)); 395 PetscCall(DMRestoreLocalVector(da, &localU)); 396 PetscCall(PetscDrawStringSetSize(draw, .2, .2)); 397 PetscCall(PetscDrawFlush(draw)); 398 PetscCall(PetscDrawSetPause(draw, pause)); 399 PetscCall(PetscDrawPause(draw)); 400 PetscFunctionReturn(PETSC_SUCCESS); 401 } 402 403 PetscErrorCode MyDestroy(void **ptr) 404 { 405 UserCtx *ctx = *(UserCtx **)ptr; 406 407 PetscFunctionBegin; 408 PetscCall(PetscDrawViewPortsDestroy(ctx->ports)); 409 PetscFunctionReturn(PETSC_SUCCESS); 410 } 411 412 /*TEST 413 414 test: 415 args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial 416 417 test: 418 suffix: 2 419 args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001 420 requires: x 421 422 TEST*/ 423