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