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 23 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 24 ./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 25 to generate InitialSolution.heat 26 27 */ 28 #include <petscdm.h> 29 #include <petscdmda.h> 30 #include <petscts.h> 31 #include <petscdraw.h> 32 33 /* 34 User-defined routines 35 */ 36 extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec),MyMonitor(TS,PetscInt,PetscReal,Vec,void*),MyDestroy(void**); 37 typedef struct {PetscReal kappa;PetscBool allencahn;PetscDrawViewPorts *ports;} UserCtx; 38 39 int main(int argc,char **argv) 40 { 41 TS ts; /* time integrator */ 42 Vec x,r; /* solution, residual vectors */ 43 PetscInt steps,Mx; 44 PetscErrorCode ierr; 45 DM da; 46 PetscReal dt; 47 UserCtx ctx; 48 PetscBool mymonitor; 49 PetscViewer viewer; 50 PetscBool flg; 51 52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 53 Initialize program 54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 55 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 56 ctx.kappa = 1.0; 57 ierr = PetscOptionsGetReal(NULL,NULL,"-kappa",&ctx.kappa,NULL);CHKERRQ(ierr); 58 ctx.allencahn = PETSC_FALSE; 59 ierr = PetscOptionsHasName(NULL,NULL,"-allen-cahn",&ctx.allencahn);CHKERRQ(ierr); 60 ierr = PetscOptionsHasName(NULL,NULL,"-mymonitor",&mymonitor);CHKERRQ(ierr); 61 62 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 63 Create distributed array (DMDA) to manage parallel grid and vectors 64 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 65 ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da);CHKERRQ(ierr); 66 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 67 ierr = DMSetUp(da);CHKERRQ(ierr); 68 ierr = DMDASetFieldName(da,0,"Heat equation: u");CHKERRQ(ierr); 69 ierr = DMDAGetInfo(da,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 70 dt = 1.0/(ctx.kappa*Mx*Mx); 71 72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 73 Extract global vectors from DMDA; then duplicate for remaining 74 vectors that are the same types 75 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 76 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr); 77 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 78 79 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 80 Create timestepping solver context 81 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 82 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 83 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 84 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 85 ierr = TSSetRHSFunction(ts,NULL,FormFunction,&ctx);CHKERRQ(ierr); 86 87 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 88 Customize nonlinear solver 89 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 90 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 91 92 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 93 Set initial conditions 94 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 95 ierr = FormInitialSolution(da,x);CHKERRQ(ierr); 96 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 97 ierr = TSSetMaxTime(ts,.02);CHKERRQ(ierr); 98 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_INTERPOLATE);CHKERRQ(ierr); 99 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 100 101 102 if (mymonitor) { 103 ctx.ports = NULL; 104 ierr = TSMonitorSet(ts,MyMonitor,&ctx,MyDestroy);CHKERRQ(ierr); 105 } 106 107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 108 Set runtime options 109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 110 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 111 112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 113 Solve nonlinear system 114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 115 ierr = TSSolve(ts,x);CHKERRQ(ierr); 116 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 117 ierr = PetscOptionsHasName(NULL,NULL,"-square_initial",&flg);CHKERRQ(ierr); 118 if (flg) { 119 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 120 ierr = VecView(x,viewer);CHKERRQ(ierr); 121 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 122 } 123 124 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 125 Free work space. All PETSc objects should be destroyed when they 126 are no longer needed. 127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 128 ierr = VecDestroy(&x);CHKERRQ(ierr); 129 ierr = VecDestroy(&r);CHKERRQ(ierr); 130 ierr = TSDestroy(&ts);CHKERRQ(ierr); 131 ierr = DMDestroy(&da);CHKERRQ(ierr); 132 133 ierr = PetscFinalize(); 134 return ierr; 135 } 136 /* ------------------------------------------------------------------- */ 137 /* 138 FormFunction - Evaluates nonlinear function, F(x). 139 140 Input Parameters: 141 . ts - the TS context 142 . X - input vector 143 . ptr - optional user-defined context, as set by SNESSetFunction() 144 145 Output Parameter: 146 . F - function vector 147 */ 148 PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr) 149 { 150 DM da; 151 PetscErrorCode ierr; 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 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 160 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 161 ierr = 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);CHKERRQ(ierr); 162 163 hx = 1.0/(PetscReal)Mx; 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 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 172 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 173 174 /* 175 Get pointers to vector data 176 */ 177 ierr = DMDAVecGetArrayRead(da,localX,&x);CHKERRQ(ierr); 178 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 179 180 /* 181 Get local grid boundaries 182 */ 183 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 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 ierr = DMDAVecRestoreArrayRead(da,localX,&x);CHKERRQ(ierr); 197 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 198 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 202 /* ------------------------------------------------------------------- */ 203 PetscErrorCode FormInitialSolution(DM da,Vec U) 204 { 205 PetscErrorCode ierr; 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 ierr = 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);CHKERRQ(ierr); 216 217 hx = 1.0/(PetscReal)Mx; 218 219 /* 220 Get pointers to vector data 221 */ 222 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 223 224 /* 225 Get local grid boundaries 226 */ 227 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 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 ierr = PetscOptionsHasName(NULL,NULL,"-square_initial",&flg);CHKERRQ(ierr); 233 if (!flg) { 234 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_READ,&viewer);CHKERRQ(ierr); 235 ierr = VecCreate(PETSC_COMM_WORLD,&finesolution);CHKERRQ(ierr); 236 ierr = VecLoad(finesolution,viewer);CHKERRQ(ierr); 237 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 238 ierr = VecGetSize(finesolution,&N);CHKERRQ(ierr); 239 scale = N/Mx; 240 ierr = VecGetArrayRead(finesolution,&f);CHKERRQ(ierr); 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 = PetscSqrtScalar((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 ierr = VecRestoreArrayRead(finesolution,&f);CHKERRQ(ierr); 260 ierr = VecDestroy(&finesolution);CHKERRQ(ierr); 261 } 262 263 /* 264 Restore vectors 265 */ 266 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 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 PetscErrorCode ierr; 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 ierr = PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,vbounds);CHKERRQ(ierr); 292 ierr = PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1200,800);CHKERRQ(ierr); 293 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 294 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 295 ierr = 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);CHKERRQ(ierr); 296 ierr = DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr); 297 hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx); 298 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 299 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 300 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 301 302 ierr = PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,&lg);CHKERRQ(ierr); 303 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 304 ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr); 305 if (!ctx->ports) { 306 ierr = PetscDrawViewPortsCreateRect(draw,1,3,&ctx->ports);CHKERRQ(ierr); 307 } 308 ports = ctx->ports; 309 ierr = PetscDrawLGGetAxis(lg,&axis);CHKERRQ(ierr); 310 ierr = PetscDrawLGReset(lg);CHKERRQ(ierr); 311 312 xx[0] = 0.0; xx[1] = 1.0; cnt = 2; 313 ierr = PetscOptionsGetRealArray(NULL,NULL,"-zoom",xx,&cnt,NULL);CHKERRQ(ierr); 314 xs = xx[0]/hx; xm = (xx[1] - xx[0])/hx; 315 316 /* 317 Plot the energies 318 */ 319 ierr = PetscDrawLGSetDimension(lg,1 + (ctx->allencahn ? 1 : 0));CHKERRQ(ierr); 320 ierr = PetscDrawLGSetColors(lg,colors+1);CHKERRQ(ierr); 321 ierr = PetscDrawViewPortsSet(ports,2);CHKERRQ(ierr); 322 x = hx*xs; 323 for (i=xs; i<xs+xm; i++) { 324 xx[0] = xx[1] = x; 325 yy[0] = PetscRealPart(.25*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx); 326 if (ctx->allencahn) yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i])); 327 ierr = PetscDrawLGAddPoint(lg,xx,yy);CHKERRQ(ierr); 328 x += hx; 329 } 330 ierr = PetscDrawGetPause(draw,&pause);CHKERRQ(ierr); 331 ierr = PetscDrawSetPause(draw,0.0);CHKERRQ(ierr); 332 ierr = PetscDrawAxisSetLabels(axis,"Energy","","");CHKERRQ(ierr); 333 ierr = PetscDrawLGSetLegend(lg,legend);CHKERRQ(ierr); 334 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 335 336 /* 337 Plot the forces 338 */ 339 ierr = PetscDrawViewPortsSet(ports,1);CHKERRQ(ierr); 340 ierr = PetscDrawLGReset(lg);CHKERRQ(ierr); 341 x = xs*hx; 342 max = 0.; 343 for (i=xs; i<xs+xm; i++) { 344 xx[0] = xx[1] = x; 345 yy[0] = PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 346 max = PetscMax(max,PetscAbs(yy[0])); 347 if (ctx->allencahn) { 348 yy[1] = PetscRealPart(u[i] - u[i]*u[i]*u[i]); 349 max = PetscMax(max,PetscAbs(yy[1])); 350 } 351 ierr = PetscDrawLGAddPoint(lg,xx,yy);CHKERRQ(ierr); 352 x += hx; 353 } 354 ierr = PetscDrawAxisSetLabels(axis,"Right hand side","","");CHKERRQ(ierr); 355 ierr = PetscDrawLGSetLegend(lg,NULL);CHKERRQ(ierr); 356 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 357 358 /* 359 Plot the solution 360 */ 361 ierr = PetscDrawLGSetDimension(lg,1);CHKERRQ(ierr); 362 ierr = PetscDrawViewPortsSet(ports,0);CHKERRQ(ierr); 363 ierr = PetscDrawLGReset(lg);CHKERRQ(ierr); 364 x = hx*xs; 365 ierr = PetscDrawLGSetLimits(lg,x,x+(xm-1)*hx,-1.1,1.1);CHKERRQ(ierr); 366 ierr = PetscDrawLGSetColors(lg,colors);CHKERRQ(ierr); 367 for (i=xs; i<xs+xm; i++) { 368 xx[0] = x; 369 yy[0] = PetscRealPart(u[i]); 370 ierr = PetscDrawLGAddPoint(lg,xx,yy);CHKERRQ(ierr); 371 x += hx; 372 } 373 ierr = PetscDrawAxisSetLabels(axis,"Solution","","");CHKERRQ(ierr); 374 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 375 376 /* 377 Print the forces as arrows on the solution 378 */ 379 x = hx*xs; 380 cnt = xm/60; 381 cnt = (!cnt) ? 1 : cnt; 382 383 for (i=xs; i<xs+xm; i += cnt) { 384 y = PetscRealPart(u[i]); 385 len = .5*PetscRealPart(ctx->kappa*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max; 386 ierr = PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_RED);CHKERRQ(ierr); 387 if (ctx->allencahn) { 388 len = .5*PetscRealPart(u[i] - u[i]*u[i]*u[i])/max; 389 ierr = PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_BLUE);CHKERRQ(ierr); 390 } 391 x += cnt*hx; 392 } 393 ierr = DMDAVecRestoreArrayRead(da,localU,&x);CHKERRQ(ierr); 394 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 395 ierr = PetscDrawStringSetSize(draw,.2,.2);CHKERRQ(ierr); 396 ierr = PetscDrawFlush(draw);CHKERRQ(ierr); 397 ierr = PetscDrawSetPause(draw,pause);CHKERRQ(ierr); 398 ierr = PetscDrawPause(draw);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 PetscErrorCode MyDestroy(void **ptr) 403 { 404 UserCtx *ctx = *(UserCtx**)ptr; 405 PetscErrorCode ierr; 406 407 PetscFunctionBegin; 408 ierr = PetscDrawViewPortsDestroy(ctx->ports);CHKERRQ(ierr); 409 PetscFunctionReturn(0); 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