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