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