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