1 2 static char help[] = "Solves heat equation in 1d.\n"; 3 4 /* 5 Solves the equation 6 7 u_t = kappa \Delta u 8 Periodic boundary conditions 9 10 Evolve the heat equation: 11 --------------- 12 ./heat -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor 13 14 Evolve the Allen-Cahn equation: 15 --------------- 16 ./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 17 18 Evolve the Allen-Cahn equation: zoom in on part of the domain 19 --------------- 20 ./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 21 22 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 23 ./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 24 to generate InitialSolution.heat 25 26 */ 27 #include <petscdm.h> 28 #include <petscdmda.h> 29 #include <petscts.h> 30 #include <petscdraw.h> 31 32 /* 33 User-defined routines 34 */ 35 extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec), MyMonitor(TS, PetscInt, PetscReal, Vec, void *), MyDestroy(void **); 36 typedef struct { 37 PetscReal kappa; 38 PetscBool allencahn; 39 PetscDrawViewPorts *ports; 40 } UserCtx; 41 42 int main(int argc, char **argv) { 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 DM da; 151 PetscInt i, Mx, xs, xm; 152 PetscReal hx, sx; 153 PetscScalar *x, *f; 154 Vec localX; 155 UserCtx *ctx = (UserCtx *)ptr; 156 157 PetscFunctionBegin; 158 PetscCall(TSGetDM(ts, &da)); 159 PetscCall(DMGetLocalVector(da, &localX)); 160 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)); 161 162 hx = 1.0 / (PetscReal)Mx; 163 sx = 1.0 / (hx * hx); 164 165 /* 166 Scatter ghost points to local vector,using the 2-step process 167 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 168 By placing code between these two statements, computations can be 169 done while messages are in transition. 170 */ 171 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 172 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 173 174 /* 175 Get pointers to vector data 176 */ 177 PetscCall(DMDAVecGetArrayRead(da, localX, &x)); 178 PetscCall(DMDAVecGetArray(da, F, &f)); 179 180 /* 181 Get local grid boundaries 182 */ 183 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 184 185 /* 186 Compute function over the locally owned part of the grid 187 */ 188 for (i = xs; i < xs + xm; i++) { 189 f[i] = ctx->kappa * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx; 190 if (ctx->allencahn) f[i] += (x[i] - x[i] * x[i] * x[i]); 191 } 192 193 /* 194 Restore vectors 195 */ 196 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x)); 197 PetscCall(DMDAVecRestoreArray(da, F, &f)); 198 PetscCall(DMRestoreLocalVector(da, &localX)); 199 PetscFunctionReturn(0); 200 } 201 202 /* ------------------------------------------------------------------- */ 203 PetscErrorCode FormInitialSolution(DM da, Vec U) { 204 PetscInt i, xs, xm, Mx, scale = 1, N; 205 PetscScalar *u; 206 const PetscScalar *f; 207 PetscReal hx, x, r; 208 Vec finesolution; 209 PetscViewer viewer; 210 PetscBool flg; 211 212 PetscFunctionBegin; 213 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)); 214 215 hx = 1.0 / (PetscReal)Mx; 216 217 /* 218 Get pointers to vector data 219 */ 220 PetscCall(DMDAVecGetArray(da, U, &u)); 221 222 /* 223 Get local grid boundaries 224 */ 225 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 226 227 /* InitialSolution is obtained with 228 ./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 229 */ 230 PetscCall(PetscOptionsHasName(NULL, NULL, "-square_initial", &flg)); 231 if (!flg) { 232 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer)); 233 PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution)); 234 PetscCall(VecLoad(finesolution, viewer)); 235 PetscCall(PetscViewerDestroy(&viewer)); 236 PetscCall(VecGetSize(finesolution, &N)); 237 scale = N / Mx; 238 PetscCall(VecGetArrayRead(finesolution, &f)); 239 } 240 241 /* 242 Compute function over the locally owned part of the grid 243 */ 244 for (i = xs; i < xs + xm; i++) { 245 x = i * hx; 246 r = PetscSqrtScalar((x - .5) * (x - .5)); 247 if (r < .125) u[i] = 1.0; 248 else u[i] = -.5; 249 250 /* With the initial condition above the method is first order in space */ 251 /* this is a smooth initial condition so the method becomes second order in space */ 252 /*u[i] = PetscSinScalar(2*PETSC_PI*x); */ 253 /* u[i] = f[scale*i];*/ 254 if (!flg) u[i] = f[scale * i]; 255 } 256 if (!flg) { 257 PetscCall(VecRestoreArrayRead(finesolution, &f)); 258 PetscCall(VecDestroy(&finesolution)); 259 } 260 261 /* 262 Restore vectors 263 */ 264 PetscCall(DMDAVecRestoreArray(da, U, &u)); 265 PetscFunctionReturn(0); 266 } 267 268 /* 269 This routine is not parallel 270 */ 271 PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr) { 272 UserCtx *ctx = (UserCtx *)ptr; 273 PetscDrawLG lg; 274 PetscScalar *u; 275 PetscInt Mx, i, xs, xm, cnt; 276 PetscReal x, y, hx, pause, sx, len, max, xx[2], yy[2]; 277 PetscDraw draw; 278 Vec localU; 279 DM da; 280 int colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE}; 281 const char *const legend[] = {"-kappa (\\grad u,\\grad u)", "(1 - u^2)^2"}; 282 PetscDrawAxis axis; 283 PetscDrawViewPorts *ports; 284 PetscReal vbounds[] = {-1.1, 1.1}; 285 286 PetscFunctionBegin; 287 PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds)); 288 PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1200, 800)); 289 PetscCall(TSGetDM(ts, &da)); 290 PetscCall(DMGetLocalVector(da, &localU)); 291 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)); 292 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 293 hx = 1.0 / (PetscReal)Mx; 294 sx = 1.0 / (hx * hx); 295 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU)); 296 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU)); 297 PetscCall(DMDAVecGetArrayRead(da, localU, &u)); 298 299 PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg)); 300 PetscCall(PetscDrawLGGetDraw(lg, &draw)); 301 PetscCall(PetscDrawCheckResizedWindow(draw)); 302 if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports)); 303 ports = ctx->ports; 304 PetscCall(PetscDrawLGGetAxis(lg, &axis)); 305 PetscCall(PetscDrawLGReset(lg)); 306 307 xx[0] = 0.0; 308 xx[1] = 1.0; 309 cnt = 2; 310 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL)); 311 xs = xx[0] / hx; 312 xm = (xx[1] - xx[0]) / hx; 313 314 /* 315 Plot the energies 316 */ 317 PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->allencahn ? 1 : 0))); 318 PetscCall(PetscDrawLGSetColors(lg, colors + 1)); 319 PetscCall(PetscDrawViewPortsSet(ports, 2)); 320 x = hx * xs; 321 for (i = xs; i < xs + xm; i++) { 322 xx[0] = xx[1] = x; 323 yy[0] = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx); 324 if (ctx->allencahn) yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i])); 325 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 326 x += hx; 327 } 328 PetscCall(PetscDrawGetPause(draw, &pause)); 329 PetscCall(PetscDrawSetPause(draw, 0.0)); 330 PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", "")); 331 PetscCall(PetscDrawLGSetLegend(lg, legend)); 332 PetscCall(PetscDrawLGDraw(lg)); 333 334 /* 335 Plot the forces 336 */ 337 PetscCall(PetscDrawViewPortsSet(ports, 1)); 338 PetscCall(PetscDrawLGReset(lg)); 339 x = xs * hx; 340 max = 0.; 341 for (i = xs; i < xs + xm; i++) { 342 xx[0] = xx[1] = x; 343 yy[0] = PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx); 344 max = PetscMax(max, PetscAbs(yy[0])); 345 if (ctx->allencahn) { 346 yy[1] = PetscRealPart(u[i] - u[i] * u[i] * u[i]); 347 max = PetscMax(max, PetscAbs(yy[1])); 348 } 349 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 350 x += hx; 351 } 352 PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", "")); 353 PetscCall(PetscDrawLGSetLegend(lg, NULL)); 354 PetscCall(PetscDrawLGDraw(lg)); 355 356 /* 357 Plot the solution 358 */ 359 PetscCall(PetscDrawLGSetDimension(lg, 1)); 360 PetscCall(PetscDrawViewPortsSet(ports, 0)); 361 PetscCall(PetscDrawLGReset(lg)); 362 x = hx * xs; 363 PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1)); 364 PetscCall(PetscDrawLGSetColors(lg, colors)); 365 for (i = xs; i < xs + xm; i++) { 366 xx[0] = x; 367 yy[0] = PetscRealPart(u[i]); 368 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 369 x += hx; 370 } 371 PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", "")); 372 PetscCall(PetscDrawLGDraw(lg)); 373 374 /* 375 Print the forces as arrows on the solution 376 */ 377 x = hx * xs; 378 cnt = xm / 60; 379 cnt = (!cnt) ? 1 : cnt; 380 381 for (i = xs; i < xs + xm; i += cnt) { 382 y = PetscRealPart(u[i]); 383 len = .5 * PetscRealPart(ctx->kappa * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max; 384 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED)); 385 if (ctx->allencahn) { 386 len = .5 * PetscRealPart(u[i] - u[i] * u[i] * u[i]) / max; 387 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE)); 388 } 389 x += cnt * hx; 390 } 391 PetscCall(DMDAVecRestoreArrayRead(da, localU, &x)); 392 PetscCall(DMRestoreLocalVector(da, &localU)); 393 PetscCall(PetscDrawStringSetSize(draw, .2, .2)); 394 PetscCall(PetscDrawFlush(draw)); 395 PetscCall(PetscDrawSetPause(draw, pause)); 396 PetscCall(PetscDrawPause(draw)); 397 PetscFunctionReturn(0); 398 } 399 400 PetscErrorCode MyDestroy(void **ptr) { 401 UserCtx *ctx = *(UserCtx **)ptr; 402 403 PetscFunctionBegin; 404 PetscCall(PetscDrawViewPortsDestroy(ctx->ports)); 405 PetscFunctionReturn(0); 406 } 407 408 /*TEST 409 410 test: 411 args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 2 -square_initial 412 413 test: 414 suffix: 2 415 args: -ts_monitor -snes_monitor -pc_type lu -snes_converged_reason -ts_type cn -da_refine 5 -mymonitor -square_initial -allen-cahn -kappa .001 416 requires: x 417 418 TEST*/ 419