1 2 static char help[] = "Solves biharmonic equation in 1d.\n"; 3 4 /* 5 Solves the equation 6 7 u_t = - kappa \Delta \Delta u 8 Periodic boundary conditions 9 10 Evolve the biharmonic heat equation: 11 --------------- 12 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor 13 14 Evolve with the restriction that -1 <= u <= 1; i.e. as a variational inequality 15 --------------- 16 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor 17 18 u_t = kappa \Delta \Delta u + 6.*u*(u_x)^2 + (3*u^2 - 12) \Delta u 19 -1 <= u <= 1 20 Periodic boundary conditions 21 22 Evolve the Cahn-Hillard equations: double well Initial hump shrinks then grows 23 --------------- 24 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -ts_monitor_draw_solution --mymonitor 25 26 Initial hump neither shrinks nor grows when degenerate (otherwise similar solution) 27 28 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -degenerate -ts_monitor_draw_solution --mymonitor 29 30 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 6 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -snes_vi_ignore_function_sign -ts_monitor_draw_solution --mymonitor 31 32 Evolve the Cahn-Hillard equations: double obstacle 33 --------------- 34 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 2 -snes_linesearch_monitor -ts_monitor_draw_solution --mymonitor 35 36 Evolve the Cahn-Hillard equations: logarithmic + double well (never shrinks and then grows) 37 --------------- 38 ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -mymonitor 39 40 ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .0001 -ts_dt 5.96046e-06 -cahn-hillard -energy 3 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --ts_max_time 1. -degenerate -mymonitor 41 42 Evolve the Cahn-Hillard equations: logarithmic + double obstacle (never shrinks, never grows) 43 --------------- 44 ./biharmonic -ts_monitor -snes_monitor -pc_type lu --snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -kappa .00001 -ts_dt 5.96046e-06 -cahn-hillard -energy 4 -snes_linesearch_monitor -theta .00000001 -ts_monitor_draw_solution --mymonitor 45 46 */ 47 #include <petscdm.h> 48 #include <petscdmda.h> 49 #include <petscts.h> 50 #include <petscdraw.h> 51 52 extern PetscErrorCode FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DM,Vec),MyMonitor(TS,PetscInt,PetscReal,Vec,void*),MyDestroy(void**),FormJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 53 typedef struct {PetscBool cahnhillard;PetscBool degenerate;PetscReal kappa;PetscInt energy;PetscReal tol;PetscReal theta,theta_c;PetscInt truncation;PetscBool netforce; PetscDrawViewPorts *ports;} UserCtx; 54 55 int main(int argc,char **argv) 56 { 57 TS ts; /* nonlinear solver */ 58 Vec x,r; /* solution, residual vectors */ 59 Mat J; /* Jacobian matrix */ 60 PetscInt steps,Mx; 61 PetscErrorCode ierr; 62 DM da; 63 PetscReal dt; 64 PetscBool mymonitor; 65 UserCtx ctx; 66 67 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 68 Initialize program 69 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 70 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 71 ctx.kappa = 1.0; 72 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-kappa",&ctx.kappa,NULL)); 73 ctx.degenerate = PETSC_FALSE; 74 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-degenerate",&ctx.degenerate,NULL)); 75 ctx.cahnhillard = PETSC_FALSE; 76 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-cahn-hillard",&ctx.cahnhillard,NULL)); 77 ctx.netforce = PETSC_FALSE; 78 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-netforce",&ctx.netforce,NULL)); 79 ctx.energy = 1; 80 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-energy",&ctx.energy,NULL)); 81 ctx.tol = 1.0e-8; 82 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-tol",&ctx.tol,NULL)); 83 ctx.theta = .001; 84 ctx.theta_c = 1.0; 85 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-theta",&ctx.theta,NULL)); 86 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-theta_c",&ctx.theta_c,NULL)); 87 ctx.truncation = 1; 88 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-truncation",&ctx.truncation,NULL)); 89 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-mymonitor",&mymonitor)); 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Create distributed array (DMDA) to manage parallel grid and vectors 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 CHKERRQ(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da)); 95 CHKERRQ(DMSetFromOptions(da)); 96 CHKERRQ(DMSetUp(da)); 97 CHKERRQ(DMDASetFieldName(da,0,"Biharmonic heat equation: u")); 98 CHKERRQ(DMDAGetInfo(da,0,&Mx,0,0,0,0,0,0,0,0,0,0,0)); 99 dt = 1.0/(10.*ctx.kappa*Mx*Mx*Mx*Mx); 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Extract global vectors from DMDA; then duplicate for remaining 103 vectors that are the same types 104 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 105 CHKERRQ(DMCreateGlobalVector(da,&x)); 106 CHKERRQ(VecDuplicate(x,&r)); 107 108 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 109 Create timestepping solver context 110 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 111 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 112 CHKERRQ(TSSetDM(ts,da)); 113 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 114 CHKERRQ(TSSetRHSFunction(ts,NULL,FormFunction,&ctx)); 115 CHKERRQ(DMSetMatType(da,MATAIJ)); 116 CHKERRQ(DMCreateMatrix(da,&J)); 117 CHKERRQ(TSSetRHSJacobian(ts,J,J,FormJacobian,&ctx)); 118 CHKERRQ(TSSetMaxTime(ts,.02)); 119 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_INTERPOLATE)); 120 121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 122 Create matrix data structure; set Jacobian evaluation routine 123 124 Set Jacobian matrix data structure and default Jacobian evaluation 125 routine. User can override with: 126 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 127 (unless user explicitly sets preconditioner) 128 -snes_mf_operator : form preconditioning matrix as set by the user, 129 but use matrix-free approx for Jacobian-vector 130 products within Newton-Krylov method 131 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 134 Customize nonlinear solver 135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 136 CHKERRQ(TSSetType(ts,TSCN)); 137 138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139 Set initial conditions 140 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 141 CHKERRQ(FormInitialSolution(da,x)); 142 CHKERRQ(TSSetTimeStep(ts,dt)); 143 CHKERRQ(TSSetSolution(ts,x)); 144 145 if (mymonitor) { 146 ctx.ports = NULL; 147 CHKERRQ(TSMonitorSet(ts,MyMonitor,&ctx,MyDestroy)); 148 } 149 150 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 151 Set runtime options 152 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 153 CHKERRQ(TSSetFromOptions(ts)); 154 155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156 Solve nonlinear system 157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 158 CHKERRQ(TSSolve(ts,x)); 159 CHKERRQ(TSGetStepNumber(ts,&steps)); 160 CHKERRQ(VecView(x,PETSC_VIEWER_BINARY_WORLD)); 161 162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 163 Free work space. All PETSc objects should be destroyed when they 164 are no longer needed. 165 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 166 CHKERRQ(MatDestroy(&J)); 167 CHKERRQ(VecDestroy(&x)); 168 CHKERRQ(VecDestroy(&r)); 169 CHKERRQ(TSDestroy(&ts)); 170 CHKERRQ(DMDestroy(&da)); 171 172 ierr = PetscFinalize(); 173 return ierr; 174 } 175 /* ------------------------------------------------------------------- */ 176 /* 177 FormFunction - Evaluates nonlinear function, F(x). 178 179 Input Parameters: 180 . ts - the TS context 181 . X - input vector 182 . ptr - optional user-defined context, as set by SNESSetFunction() 183 184 Output Parameter: 185 . F - function vector 186 */ 187 PetscErrorCode FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void *ptr) 188 { 189 DM da; 190 PetscInt i,Mx,xs,xm; 191 PetscReal hx,sx; 192 PetscScalar *x,*f,c,r,l; 193 Vec localX; 194 UserCtx *ctx = (UserCtx*)ptr; 195 PetscReal tol = ctx->tol, theta=ctx->theta,theta_c=ctx->theta_c,a,b; /* a and b are used in the cubic truncation of the log function */ 196 197 PetscFunctionBegin; 198 CHKERRQ(TSGetDM(ts,&da)); 199 CHKERRQ(DMGetLocalVector(da,&localX)); 200 CHKERRQ(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)); 201 202 hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx); 203 204 /* 205 Scatter ghost points to local vector,using the 2-step process 206 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 207 By placing code between these two statements, computations can be 208 done while messages are in transition. 209 */ 210 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 211 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 212 213 /* 214 Get pointers to vector data 215 */ 216 CHKERRQ(DMDAVecGetArrayRead(da,localX,&x)); 217 CHKERRQ(DMDAVecGetArray(da,F,&f)); 218 219 /* 220 Get local grid boundaries 221 */ 222 CHKERRQ(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 223 224 /* 225 Compute function over the locally owned part of the grid 226 */ 227 for (i=xs; i<xs+xm; i++) { 228 if (ctx->degenerate) { 229 c = (1. - x[i]*x[i])*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 230 r = (1. - x[i+1]*x[i+1])*(x[i] + x[i+2] - 2.0*x[i+1])*sx; 231 l = (1. - x[i-1]*x[i-1])*(x[i-2] + x[i] - 2.0*x[i-1])*sx; 232 } else { 233 c = (x[i-1] + x[i+1] - 2.0*x[i])*sx; 234 r = (x[i] + x[i+2] - 2.0*x[i+1])*sx; 235 l = (x[i-2] + x[i] - 2.0*x[i-1])*sx; 236 } 237 f[i] = -ctx->kappa*(l + r - 2.0*c)*sx; 238 if (ctx->cahnhillard) { 239 switch (ctx->energy) { 240 case 1: /* double well */ 241 f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 242 break; 243 case 2: /* double obstacle */ 244 f[i] += -(x[i-1] + x[i+1] - 2.0*x[i])*sx; 245 break; 246 case 3: /* logarithmic + double well */ 247 f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 248 if (ctx->truncation==2) { /* log function with approximated with a quadratic polynomial outside -1.0+2*tol, 1.0-2*tol */ 249 if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 250 else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 251 else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 252 } else { /* log function is approximated with a cubic polynomial outside -1.0+2*tol, 1.0-2*tol */ 253 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 254 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 255 if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += -1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (-1.0*a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 256 else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += 1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + ( a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 257 else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 258 } 259 break; 260 case 4: /* logarithmic + double obstacle */ 261 f[i] += -theta_c*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 262 if (ctx->truncation==2) { /* quadratic */ 263 if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 264 else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += (.25*theta/(tol-tol*tol))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 265 else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 266 } else { /* cubic */ 267 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 268 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 269 if (PetscRealPart(x[i]) < -1.0 + 2.0*tol) f[i] += -1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (-1.0*a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 270 else if (PetscRealPart(x[i]) > 1.0 - 2.0*tol) f[i] += 1.0*a*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + ( a*x[i] + b)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 271 else f[i] += 2.0*theta*x[i]/((1.0-x[i]*x[i])*(1.0-x[i]*x[i]))*.25*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (theta/(1.0-x[i]*x[i]))*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 272 } 273 break; 274 } 275 } 276 277 } 278 279 /* 280 Restore vectors 281 */ 282 CHKERRQ(DMDAVecRestoreArrayRead(da,localX,&x)); 283 CHKERRQ(DMDAVecRestoreArray(da,F,&f)); 284 CHKERRQ(DMRestoreLocalVector(da,&localX)); 285 PetscFunctionReturn(0); 286 } 287 288 /* ------------------------------------------------------------------- */ 289 /* 290 FormJacobian - Evaluates nonlinear function's Jacobian 291 292 */ 293 PetscErrorCode FormJacobian(TS ts,PetscReal ftime,Vec X,Mat A,Mat B,void *ptr) 294 { 295 DM da; 296 PetscInt i,Mx,xs,xm; 297 MatStencil row,cols[5]; 298 PetscReal hx,sx; 299 PetscScalar *x,vals[5]; 300 Vec localX; 301 UserCtx *ctx = (UserCtx*)ptr; 302 303 PetscFunctionBegin; 304 CHKERRQ(TSGetDM(ts,&da)); 305 CHKERRQ(DMGetLocalVector(da,&localX)); 306 CHKERRQ(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)); 307 308 hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx); 309 310 /* 311 Scatter ghost points to local vector,using the 2-step process 312 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 313 By placing code between these two statements, computations can be 314 done while messages are in transition. 315 */ 316 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 317 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 318 319 /* 320 Get pointers to vector data 321 */ 322 CHKERRQ(DMDAVecGetArrayRead(da,localX,&x)); 323 324 /* 325 Get local grid boundaries 326 */ 327 CHKERRQ(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 328 329 /* 330 Compute function over the locally owned part of the grid 331 */ 332 for (i=xs; i<xs+xm; i++) { 333 row.i = i; 334 if (ctx->degenerate) { 335 /*PetscScalar c,r,l; 336 c = (1. - x[i]*x[i])*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 337 r = (1. - x[i+1]*x[i+1])*(x[i] + x[i+2] - 2.0*x[i+1])*sx; 338 l = (1. - x[i-1]*x[i-1])*(x[i-2] + x[i] - 2.0*x[i-1])*sx; */ 339 } else { 340 cols[0].i = i - 2; vals[0] = -ctx->kappa*sx*sx; 341 cols[1].i = i - 1; vals[1] = 4.0*ctx->kappa*sx*sx; 342 cols[2].i = i ; vals[2] = -6.0*ctx->kappa*sx*sx; 343 cols[3].i = i + 1; vals[3] = 4.0*ctx->kappa*sx*sx; 344 cols[4].i = i + 2; vals[4] = -ctx->kappa*sx*sx; 345 } 346 CHKERRQ(MatSetValuesStencil(B,1,&row,5,cols,vals,INSERT_VALUES)); 347 348 if (ctx->cahnhillard) { 349 switch (ctx->energy) { 350 case 1: /* double well */ 351 /* f[i] += 6.*.25*x[i]*(x[i+1] - x[i-1])*(x[i+1] - x[i-1])*sx + (3.*x[i]*x[i] - 1.)*(x[i-1] + x[i+1] - 2.0*x[i])*sx; */ 352 break; 353 case 2: /* double obstacle */ 354 /* f[i] += -(x[i-1] + x[i+1] - 2.0*x[i])*sx; */ 355 break; 356 case 3: /* logarithmic + double well */ 357 break; 358 case 4: /* logarithmic + double obstacle */ 359 break; 360 } 361 } 362 363 } 364 365 /* 366 Restore vectors 367 */ 368 CHKERRQ(DMDAVecRestoreArrayRead(da,localX,&x)); 369 CHKERRQ(DMRestoreLocalVector(da,&localX)); 370 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 371 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 372 if (A != B) { 373 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 374 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 375 } 376 PetscFunctionReturn(0); 377 } 378 /* ------------------------------------------------------------------- */ 379 PetscErrorCode FormInitialSolution(DM da,Vec U) 380 { 381 PetscInt i,xs,xm,Mx,N,scale; 382 PetscScalar *u; 383 PetscReal r,hx,x; 384 const PetscScalar *f; 385 Vec finesolution; 386 PetscViewer viewer; 387 388 PetscFunctionBegin; 389 CHKERRQ(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)); 390 391 hx = 1.0/(PetscReal)Mx; 392 393 /* 394 Get pointers to vector data 395 */ 396 CHKERRQ(DMDAVecGetArray(da,U,&u)); 397 398 /* 399 Get local grid boundaries 400 */ 401 CHKERRQ(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 402 403 /* 404 Seee heat.c for how to generate InitialSolution.heat 405 */ 406 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_READ,&viewer)); 407 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&finesolution)); 408 CHKERRQ(VecLoad(finesolution,viewer)); 409 CHKERRQ(PetscViewerDestroy(&viewer)); 410 CHKERRQ(VecGetSize(finesolution,&N)); 411 scale = N/Mx; 412 CHKERRQ(VecGetArrayRead(finesolution,&f)); 413 414 /* 415 Compute function over the locally owned part of the grid 416 */ 417 for (i=xs; i<xs+xm; i++) { 418 x = i*hx; 419 r = PetscSqrtReal((x-.5)*(x-.5)); 420 if (r < .125) u[i] = 1.0; 421 else u[i] = -.5; 422 423 /* With the initial condition above the method is first order in space */ 424 /* this is a smooth initial condition so the method becomes second order in space */ 425 /*u[i] = PetscSinScalar(2*PETSC_PI*x); */ 426 u[i] = f[scale*i]; 427 } 428 CHKERRQ(VecRestoreArrayRead(finesolution,&f)); 429 CHKERRQ(VecDestroy(&finesolution)); 430 431 /* 432 Restore vectors 433 */ 434 CHKERRQ(DMDAVecRestoreArray(da,U,&u)); 435 PetscFunctionReturn(0); 436 } 437 438 /* 439 This routine is not parallel 440 */ 441 PetscErrorCode MyMonitor(TS ts,PetscInt step,PetscReal time,Vec U,void *ptr) 442 { 443 UserCtx *ctx = (UserCtx*)ptr; 444 PetscDrawLG lg; 445 PetscErrorCode ierr; 446 PetscScalar *u,l,r,c; 447 PetscInt Mx,i,xs,xm,cnt; 448 PetscReal x,y,hx,pause,sx,len,max,xx[4],yy[4],xx_netforce,yy_netforce,yup,ydown,y2,len2; 449 PetscDraw draw; 450 Vec localU; 451 DM da; 452 int colors[] = {PETSC_DRAW_YELLOW,PETSC_DRAW_RED,PETSC_DRAW_BLUE,PETSC_DRAW_PLUM,PETSC_DRAW_BLACK}; 453 /* 454 const char *const legend[3][3] = {{"-kappa (\\grad u,\\grad u)","(1 - u^2)^2"},{"-kappa (\\grad u,\\grad u)","(1 - u^2)"},{"-kappa (\\grad u,\\grad u)","logarithmic"}}; 455 */ 456 PetscDrawAxis axis; 457 PetscDrawViewPorts *ports; 458 PetscReal tol = ctx->tol, theta=ctx->theta,theta_c=ctx->theta_c,a,b; /* a and b are used in the cubic truncation of the log function */ 459 PetscReal vbounds[] = {-1.1,1.1}; 460 461 PetscFunctionBegin; 462 CHKERRQ(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,vbounds)); 463 CHKERRQ(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),800,600)); 464 CHKERRQ(TSGetDM(ts,&da)); 465 CHKERRQ(DMGetLocalVector(da,&localU)); 466 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 467 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 468 CHKERRQ(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 469 hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx); 470 CHKERRQ(DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU)); 471 CHKERRQ(DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU)); 472 CHKERRQ(DMDAVecGetArrayRead(da,localU,&u)); 473 474 CHKERRQ(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,&lg)); 475 CHKERRQ(PetscDrawLGGetDraw(lg,&draw)); 476 CHKERRQ(PetscDrawCheckResizedWindow(draw)); 477 if (!ctx->ports) { 478 CHKERRQ(PetscDrawViewPortsCreateRect(draw,1,3,&ctx->ports)); 479 } 480 ports = ctx->ports; 481 CHKERRQ(PetscDrawLGGetAxis(lg,&axis)); 482 CHKERRQ(PetscDrawLGReset(lg)); 483 484 xx[0] = 0.0; xx[1] = 1.0; cnt = 2; 485 CHKERRQ(PetscOptionsGetRealArray(NULL,NULL,"-zoom",xx,&cnt,NULL)); 486 xs = xx[0]/hx; xm = (xx[1] - xx[0])/hx; 487 488 /* 489 Plot the energies 490 */ 491 CHKERRQ(PetscDrawLGSetDimension(lg,1 + (ctx->cahnhillard ? 1 : 0) + (ctx->energy == 3))); 492 CHKERRQ(PetscDrawLGSetColors(lg,colors+1)); 493 CHKERRQ(PetscDrawViewPortsSet(ports,2)); 494 x = hx*xs; 495 for (i=xs; i<xs+xm; i++) { 496 xx[0] = xx[1] = xx[2] = x; 497 if (ctx->degenerate) yy[0] = PetscRealPart(.25*(1. - u[i]*u[i])*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx); 498 else yy[0] = PetscRealPart(.25*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx); 499 500 if (ctx->cahnhillard) { 501 switch (ctx->energy) { 502 case 1: /* double well */ 503 yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i])); 504 break; 505 case 2: /* double obstacle */ 506 yy[1] = .5*PetscRealPart(1. - u[i]*u[i]); 507 break; 508 case 3: /* logarithm + double well */ 509 yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i])); 510 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = .5*theta*(2.0*tol*PetscLogReal(tol) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.-u[i])/2.0)); 511 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + 2.0*tol*PetscLogReal(tol)); 512 else yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.0-u[i])/2.0)); 513 break; 514 case 4: /* logarithm + double obstacle */ 515 yy[1] = .5*theta_c*PetscRealPart(1.0-u[i]*u[i]); 516 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = .5*theta*(2.0*tol*PetscLogReal(tol) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.-u[i])/2.0)); 517 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + 2.0*tol*PetscLogReal(tol)); 518 else yy[2] = .5*theta*(PetscRealPart(1.0+u[i])*PetscLogReal(PetscRealPart(1.0+u[i])/2.0) + PetscRealPart(1.0-u[i])*PetscLogReal(PetscRealPart(1.0-u[i])/2.0)); 519 break; 520 default: 521 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"It will always be one of the values"); 522 } 523 } 524 CHKERRQ(PetscDrawLGAddPoint(lg,xx,yy)); 525 x += hx; 526 } 527 CHKERRQ(PetscDrawGetPause(draw,&pause)); 528 CHKERRQ(PetscDrawSetPause(draw,0.0)); 529 CHKERRQ(PetscDrawAxisSetLabels(axis,"Energy","","")); 530 /* CHKERRQ(PetscDrawLGSetLegend(lg,legend[ctx->energy-1])); */ 531 CHKERRQ(PetscDrawLGDraw(lg)); 532 533 /* 534 Plot the forces 535 */ 536 CHKERRQ(PetscDrawLGSetDimension(lg,0 + (ctx->cahnhillard ? 2 : 0) + (ctx->energy == 3))); 537 CHKERRQ(PetscDrawLGSetColors(lg,colors+1)); 538 CHKERRQ(PetscDrawViewPortsSet(ports,1)); 539 CHKERRQ(PetscDrawLGReset(lg)); 540 x = xs*hx; 541 max = 0.; 542 for (i=xs; i<xs+xm; i++) { 543 xx[0] = xx[1] = xx[2] = xx[3] = x; 544 xx_netforce = x; 545 if (ctx->degenerate) { 546 c = (1. - u[i]*u[i])*(u[i-1] + u[i+1] - 2.0*u[i])*sx; 547 r = (1. - u[i+1]*u[i+1])*(u[i] + u[i+2] - 2.0*u[i+1])*sx; 548 l = (1. - u[i-1]*u[i-1])*(u[i-2] + u[i] - 2.0*u[i-1])*sx; 549 } else { 550 c = (u[i-1] + u[i+1] - 2.0*u[i])*sx; 551 r = (u[i] + u[i+2] - 2.0*u[i+1])*sx; 552 l = (u[i-2] + u[i] - 2.0*u[i-1])*sx; 553 } 554 yy[0] = PetscRealPart(-ctx->kappa*(l + r - 2.0*c)*sx); 555 yy_netforce = yy[0]; 556 max = PetscMax(max,PetscAbs(yy[0])); 557 if (ctx->cahnhillard) { 558 switch (ctx->energy) { 559 case 1: /* double well */ 560 yy[1] = PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 561 break; 562 case 2: /* double obstacle */ 563 yy[1] = -PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx; 564 break; 565 case 3: /* logarithmic + double well */ 566 yy[1] = PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 567 if (ctx->truncation==2) { /* quadratic */ 568 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx; 569 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx; 570 else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 571 } else { /* cubic */ 572 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 573 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 574 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = PetscRealPart(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 575 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = PetscRealPart(1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 576 else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 577 } 578 break; 579 case 4: /* logarithmic + double obstacle */ 580 yy[1] = theta_c*PetscRealPart(-(u[i-1] + u[i+1] - 2.0*u[i]))*sx; 581 if (ctx->truncation==2) { 582 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx; 583 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = (.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx; 584 else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 585 } 586 else { 587 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 588 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 589 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) yy[2] = PetscRealPart(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 590 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) yy[2] = PetscRealPart(1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 591 else yy[2] = PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx); 592 } 593 break; 594 default: 595 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"It will always be one of the values"); 596 } 597 if (ctx->energy < 3) { 598 max = PetscMax(max,PetscAbs(yy[1])); 599 yy[2] = yy[0]+yy[1]; 600 yy_netforce = yy[2]; 601 } else { 602 max = PetscMax(max,PetscAbs(yy[1]+yy[2])); 603 yy[3] = yy[0]+yy[1]+yy[2]; 604 yy_netforce = yy[3]; 605 } 606 } 607 if (ctx->netforce) { 608 CHKERRQ(PetscDrawLGAddPoint(lg,&xx_netforce,&yy_netforce)); 609 } else { 610 CHKERRQ(PetscDrawLGAddPoint(lg,xx,yy)); 611 } 612 x += hx; 613 /*if (max > 7200150000.0) */ 614 /* printf("max very big when i = %d\n",i); */ 615 } 616 CHKERRQ(PetscDrawAxisSetLabels(axis,"Right hand side","","")); 617 CHKERRQ(PetscDrawLGSetLegend(lg,NULL)); 618 CHKERRQ(PetscDrawLGDraw(lg)); 619 620 /* 621 Plot the solution 622 */ 623 CHKERRQ(PetscDrawLGSetDimension(lg,1)); 624 CHKERRQ(PetscDrawViewPortsSet(ports,0)); 625 CHKERRQ(PetscDrawLGReset(lg)); 626 x = hx*xs; 627 CHKERRQ(PetscDrawLGSetLimits(lg,x,x+(xm-1)*hx,-1.1,1.1)); 628 CHKERRQ(PetscDrawLGSetColors(lg,colors)); 629 for (i=xs; i<xs+xm; i++) { 630 xx[0] = x; 631 yy[0] = PetscRealPart(u[i]); 632 CHKERRQ(PetscDrawLGAddPoint(lg,xx,yy)); 633 x += hx; 634 } 635 CHKERRQ(PetscDrawAxisSetLabels(axis,"Solution","","")); 636 CHKERRQ(PetscDrawLGDraw(lg)); 637 638 /* 639 Print the forces as arrows on the solution 640 */ 641 x = hx*xs; 642 cnt = xm/60; 643 cnt = (!cnt) ? 1 : cnt; 644 645 for (i=xs; i<xs+xm; i += cnt) { 646 y = yup = ydown = PetscRealPart(u[i]); 647 c = (u[i-1] + u[i+1] - 2.0*u[i])*sx; 648 r = (u[i] + u[i+2] - 2.0*u[i+1])*sx; 649 l = (u[i-2] + u[i] - 2.0*u[i-1])*sx; 650 len = -.5*PetscRealPart(ctx->kappa*(l + r - 2.0*c)*sx)/max; 651 CHKERRQ(PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_RED)); 652 if (ctx->cahnhillard) { 653 if (len < 0.) ydown += len; 654 else yup += len; 655 656 switch (ctx->energy) { 657 case 1: /* double well */ 658 len = .5*PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max; 659 break; 660 case 2: /* double obstacle */ 661 len = -.5*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max; 662 break; 663 case 3: /* logarithmic + double well */ 664 len = .5*PetscRealPart(6.*.25*u[i]*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (3.*u[i]*u[i] - 1.)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max; 665 if (len < 0.) ydown += len; 666 else yup += len; 667 668 if (ctx->truncation==2) { /* quadratic */ 669 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max; 670 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max; 671 else len2 = PetscRealPart(.5*(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max); 672 } else { /* cubic */ 673 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 674 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 675 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = PetscRealPart(.5*(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max); 676 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = PetscRealPart(.5*(a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max); 677 else len2 = PetscRealPart(.5*(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max); 678 } 679 y2 = len < 0 ? ydown : yup; 680 CHKERRQ(PetscDrawArrow(draw,x,y2,x,y2+len2,PETSC_DRAW_PLUM)); 681 break; 682 case 4: /* logarithmic + double obstacle */ 683 len = -.5*theta_c*PetscRealPart(-(u[i-1] + u[i+1] - 2.0*u[i])*sx/max); 684 if (len < 0.) ydown += len; 685 else yup += len; 686 687 if (ctx->truncation==2) { /* quadratic */ 688 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max; 689 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = .5*(.25*theta/(tol-tol*tol))*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max; 690 else len2 = PetscRealPart(.5*(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max); 691 } else { /* cubic */ 692 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 693 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 694 if (PetscRealPart(u[i]) < -1.0 + 2.0*tol) len2 = .5*PetscRealPart(-1.0*a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (-1.0*a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max; 695 else if (PetscRealPart(u[i]) > 1.0 - 2.0*tol) len2 = .5*PetscRealPart(a*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + ( a*u[i] + b)*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max; 696 else len2 = .5*PetscRealPart(2.0*theta*u[i]/((1.0-u[i]*u[i])*(1.0-u[i]*u[i]))*.25*(u[i+1] - u[i-1])*(u[i+1] - u[i-1])*sx + (theta/(1.0-u[i]*u[i]))*(u[i-1] + u[i+1] - 2.0*u[i])*sx)/max; 697 } 698 y2 = len < 0 ? ydown : yup; 699 CHKERRQ(PetscDrawArrow(draw,x,y2,x,y2+len2,PETSC_DRAW_PLUM)); 700 break; 701 } 702 CHKERRQ(PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_BLUE)); 703 } 704 x += cnt*hx; 705 } 706 CHKERRQ(DMDAVecRestoreArrayRead(da,localU,&x)); 707 CHKERRQ(DMRestoreLocalVector(da,&localU)); 708 CHKERRQ(PetscDrawStringSetSize(draw,.2,.2)); 709 CHKERRQ(PetscDrawFlush(draw)); 710 CHKERRQ(PetscDrawSetPause(draw,pause)); 711 CHKERRQ(PetscDrawPause(draw)); 712 PetscFunctionReturn(0); 713 } 714 715 PetscErrorCode MyDestroy(void **ptr) 716 { 717 UserCtx *ctx = *(UserCtx**)ptr; 718 719 PetscFunctionBegin; 720 CHKERRQ(PetscDrawViewPortsDestroy(ctx->ports)); 721 PetscFunctionReturn(0); 722 } 723 724 /*TEST 725 726 test: 727 TODO: currently requires initial condition file generated by heat 728 729 TEST*/ 730