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 DM da; 62 PetscReal dt; 63 PetscBool mymonitor; 64 UserCtx ctx; 65 66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 67 Initialize program 68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 69 PetscFunctionBeginUser; 70 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 71 ctx.kappa = 1.0; 72 PetscCall(PetscOptionsGetReal(NULL,NULL,"-kappa",&ctx.kappa,NULL)); 73 ctx.degenerate = PETSC_FALSE; 74 PetscCall(PetscOptionsGetBool(NULL,NULL,"-degenerate",&ctx.degenerate,NULL)); 75 ctx.cahnhillard = PETSC_FALSE; 76 PetscCall(PetscOptionsGetBool(NULL,NULL,"-cahn-hillard",&ctx.cahnhillard,NULL)); 77 ctx.netforce = PETSC_FALSE; 78 PetscCall(PetscOptionsGetBool(NULL,NULL,"-netforce",&ctx.netforce,NULL)); 79 ctx.energy = 1; 80 PetscCall(PetscOptionsGetInt(NULL,NULL,"-energy",&ctx.energy,NULL)); 81 ctx.tol = 1.0e-8; 82 PetscCall(PetscOptionsGetReal(NULL,NULL,"-tol",&ctx.tol,NULL)); 83 ctx.theta = .001; 84 ctx.theta_c = 1.0; 85 PetscCall(PetscOptionsGetReal(NULL,NULL,"-theta",&ctx.theta,NULL)); 86 PetscCall(PetscOptionsGetReal(NULL,NULL,"-theta_c",&ctx.theta_c,NULL)); 87 ctx.truncation = 1; 88 PetscCall(PetscOptionsGetInt(NULL,NULL,"-truncation",&ctx.truncation,NULL)); 89 PetscCall(PetscOptionsHasName(NULL,NULL,"-mymonitor",&mymonitor)); 90 91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 92 Create distributed array (DMDA) to manage parallel grid and vectors 93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 94 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da)); 95 PetscCall(DMSetFromOptions(da)); 96 PetscCall(DMSetUp(da)); 97 PetscCall(DMDASetFieldName(da,0,"Biharmonic heat equation: u")); 98 PetscCall(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 PetscCall(DMCreateGlobalVector(da,&x)); 106 PetscCall(VecDuplicate(x,&r)); 107 108 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 109 Create timestepping solver context 110 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 111 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 112 PetscCall(TSSetDM(ts,da)); 113 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 114 PetscCall(TSSetRHSFunction(ts,NULL,FormFunction,&ctx)); 115 PetscCall(DMSetMatType(da,MATAIJ)); 116 PetscCall(DMCreateMatrix(da,&J)); 117 PetscCall(TSSetRHSJacobian(ts,J,J,FormJacobian,&ctx)); 118 PetscCall(TSSetMaxTime(ts,.02)); 119 PetscCall(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 PetscCall(TSSetType(ts,TSCN)); 137 138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139 Set initial conditions 140 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 141 PetscCall(FormInitialSolution(da,x)); 142 PetscCall(TSSetTimeStep(ts,dt)); 143 PetscCall(TSSetSolution(ts,x)); 144 145 if (mymonitor) { 146 ctx.ports = NULL; 147 PetscCall(TSMonitorSet(ts,MyMonitor,&ctx,MyDestroy)); 148 } 149 150 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 151 Set runtime options 152 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 153 PetscCall(TSSetFromOptions(ts)); 154 155 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 156 Solve nonlinear system 157 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 158 PetscCall(TSSolve(ts,x)); 159 PetscCall(TSGetStepNumber(ts,&steps)); 160 PetscCall(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 PetscCall(MatDestroy(&J)); 167 PetscCall(VecDestroy(&x)); 168 PetscCall(VecDestroy(&r)); 169 PetscCall(TSDestroy(&ts)); 170 PetscCall(DMDestroy(&da)); 171 172 PetscCall(PetscFinalize()); 173 return 0; 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 PetscCall(TSGetDM(ts,&da)); 199 PetscCall(DMGetLocalVector(da,&localX)); 200 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)); 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 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 211 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 212 213 /* 214 Get pointers to vector data 215 */ 216 PetscCall(DMDAVecGetArrayRead(da,localX,&x)); 217 PetscCall(DMDAVecGetArray(da,F,&f)); 218 219 /* 220 Get local grid boundaries 221 */ 222 PetscCall(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 PetscCall(DMDAVecRestoreArrayRead(da,localX,&x)); 283 PetscCall(DMDAVecRestoreArray(da,F,&f)); 284 PetscCall(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 PetscCall(TSGetDM(ts,&da)); 305 PetscCall(DMGetLocalVector(da,&localX)); 306 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)); 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 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 317 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 318 319 /* 320 Get pointers to vector data 321 */ 322 PetscCall(DMDAVecGetArrayRead(da,localX,&x)); 323 324 /* 325 Get local grid boundaries 326 */ 327 PetscCall(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 PetscCall(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 PetscCall(DMDAVecRestoreArrayRead(da,localX,&x)); 369 PetscCall(DMRestoreLocalVector(da,&localX)); 370 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 371 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 372 if (A != B) { 373 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 374 PetscCall(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 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)); 390 391 hx = 1.0/(PetscReal)Mx; 392 393 /* 394 Get pointers to vector data 395 */ 396 PetscCall(DMDAVecGetArray(da,U,&u)); 397 398 /* 399 Get local grid boundaries 400 */ 401 PetscCall(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 402 403 /* 404 Seee heat.c for how to generate InitialSolution.heat 405 */ 406 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,"InitialSolution.heat",FILE_MODE_READ,&viewer)); 407 PetscCall(VecCreate(PETSC_COMM_WORLD,&finesolution)); 408 PetscCall(VecLoad(finesolution,viewer)); 409 PetscCall(PetscViewerDestroy(&viewer)); 410 PetscCall(VecGetSize(finesolution,&N)); 411 scale = N/Mx; 412 PetscCall(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 PetscCall(VecRestoreArrayRead(finesolution,&f)); 429 PetscCall(VecDestroy(&finesolution)); 430 431 /* 432 Restore vectors 433 */ 434 PetscCall(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 PetscScalar *u,l,r,c; 446 PetscInt Mx,i,xs,xm,cnt; 447 PetscReal x,y,hx,pause,sx,len,max,xx[4],yy[4],xx_netforce,yy_netforce,yup,ydown,y2,len2; 448 PetscDraw draw; 449 Vec localU; 450 DM da; 451 int colors[] = {PETSC_DRAW_YELLOW,PETSC_DRAW_RED,PETSC_DRAW_BLUE,PETSC_DRAW_PLUM,PETSC_DRAW_BLACK}; 452 /* 453 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"}}; 454 */ 455 PetscDrawAxis axis; 456 PetscDrawViewPorts *ports; 457 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 */ 458 PetscReal vbounds[] = {-1.1,1.1}; 459 460 PetscFunctionBegin; 461 PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,vbounds)); 462 PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),800,600)); 463 PetscCall(TSGetDM(ts,&da)); 464 PetscCall(DMGetLocalVector(da,&localU)); 465 PetscCall(DMDAGetInfo(da,PETSC_IGNORE,&Mx,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 466 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE)); 467 PetscCall(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL)); 468 hx = 1.0/(PetscReal)Mx; sx = 1.0/(hx*hx); 469 PetscCall(DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU)); 470 PetscCall(DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU)); 471 PetscCall(DMDAVecGetArrayRead(da,localU,&u)); 472 473 PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD),1,&lg)); 474 PetscCall(PetscDrawLGGetDraw(lg,&draw)); 475 PetscCall(PetscDrawCheckResizedWindow(draw)); 476 if (!ctx->ports) { 477 PetscCall(PetscDrawViewPortsCreateRect(draw,1,3,&ctx->ports)); 478 } 479 ports = ctx->ports; 480 PetscCall(PetscDrawLGGetAxis(lg,&axis)); 481 PetscCall(PetscDrawLGReset(lg)); 482 483 xx[0] = 0.0; xx[1] = 1.0; cnt = 2; 484 PetscCall(PetscOptionsGetRealArray(NULL,NULL,"-zoom",xx,&cnt,NULL)); 485 xs = xx[0]/hx; xm = (xx[1] - xx[0])/hx; 486 487 /* 488 Plot the energies 489 */ 490 PetscCall(PetscDrawLGSetDimension(lg,1 + (ctx->cahnhillard ? 1 : 0) + (ctx->energy == 3))); 491 PetscCall(PetscDrawLGSetColors(lg,colors+1)); 492 PetscCall(PetscDrawViewPortsSet(ports,2)); 493 x = hx*xs; 494 for (i=xs; i<xs+xm; i++) { 495 xx[0] = xx[1] = xx[2] = x; 496 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); 497 else yy[0] = PetscRealPart(.25*ctx->kappa*(u[i-1] - u[i+1])*(u[i-1] - u[i+1])*sx); 498 499 if (ctx->cahnhillard) { 500 switch (ctx->energy) { 501 case 1: /* double well */ 502 yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i])); 503 break; 504 case 2: /* double obstacle */ 505 yy[1] = .5*PetscRealPart(1. - u[i]*u[i]); 506 break; 507 case 3: /* logarithm + double well */ 508 yy[1] = .25*PetscRealPart((1. - u[i]*u[i])*(1. - u[i]*u[i])); 509 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)); 510 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)); 511 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)); 512 break; 513 case 4: /* logarithm + double obstacle */ 514 yy[1] = .5*theta_c*PetscRealPart(1.0-u[i]*u[i]); 515 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)); 516 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)); 517 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)); 518 break; 519 default: 520 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"It will always be one of the values"); 521 } 522 } 523 PetscCall(PetscDrawLGAddPoint(lg,xx,yy)); 524 x += hx; 525 } 526 PetscCall(PetscDrawGetPause(draw,&pause)); 527 PetscCall(PetscDrawSetPause(draw,0.0)); 528 PetscCall(PetscDrawAxisSetLabels(axis,"Energy","","")); 529 /* PetscCall(PetscDrawLGSetLegend(lg,legend[ctx->energy-1])); */ 530 PetscCall(PetscDrawLGDraw(lg)); 531 532 /* 533 Plot the forces 534 */ 535 PetscCall(PetscDrawLGSetDimension(lg,0 + (ctx->cahnhillard ? 2 : 0) + (ctx->energy == 3))); 536 PetscCall(PetscDrawLGSetColors(lg,colors+1)); 537 PetscCall(PetscDrawViewPortsSet(ports,1)); 538 PetscCall(PetscDrawLGReset(lg)); 539 x = xs*hx; 540 max = 0.; 541 for (i=xs; i<xs+xm; i++) { 542 xx[0] = xx[1] = xx[2] = xx[3] = x; 543 xx_netforce = x; 544 if (ctx->degenerate) { 545 c = (1. - u[i]*u[i])*(u[i-1] + u[i+1] - 2.0*u[i])*sx; 546 r = (1. - u[i+1]*u[i+1])*(u[i] + u[i+2] - 2.0*u[i+1])*sx; 547 l = (1. - u[i-1]*u[i-1])*(u[i-2] + u[i] - 2.0*u[i-1])*sx; 548 } else { 549 c = (u[i-1] + u[i+1] - 2.0*u[i])*sx; 550 r = (u[i] + u[i+2] - 2.0*u[i+1])*sx; 551 l = (u[i-2] + u[i] - 2.0*u[i-1])*sx; 552 } 553 yy[0] = PetscRealPart(-ctx->kappa*(l + r - 2.0*c)*sx); 554 yy_netforce = yy[0]; 555 max = PetscMax(max,PetscAbs(yy[0])); 556 if (ctx->cahnhillard) { 557 switch (ctx->energy) { 558 case 1: /* double well */ 559 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); 560 break; 561 case 2: /* double obstacle */ 562 yy[1] = -PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx; 563 break; 564 case 3: /* logarithmic + double well */ 565 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); 566 if (ctx->truncation==2) { /* quadratic */ 567 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; 568 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; 569 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); 570 } else { /* cubic */ 571 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 572 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 573 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); 574 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); 575 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); 576 } 577 break; 578 case 4: /* logarithmic + double obstacle */ 579 yy[1] = theta_c*PetscRealPart(-(u[i-1] + u[i+1] - 2.0*u[i]))*sx; 580 if (ctx->truncation==2) { 581 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; 582 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; 583 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); 584 } 585 else { 586 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 587 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 588 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); 589 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); 590 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); 591 } 592 break; 593 default: 594 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"It will always be one of the values"); 595 } 596 if (ctx->energy < 3) { 597 max = PetscMax(max,PetscAbs(yy[1])); 598 yy[2] = yy[0]+yy[1]; 599 yy_netforce = yy[2]; 600 } else { 601 max = PetscMax(max,PetscAbs(yy[1]+yy[2])); 602 yy[3] = yy[0]+yy[1]+yy[2]; 603 yy_netforce = yy[3]; 604 } 605 } 606 if (ctx->netforce) { 607 PetscCall(PetscDrawLGAddPoint(lg,&xx_netforce,&yy_netforce)); 608 } else { 609 PetscCall(PetscDrawLGAddPoint(lg,xx,yy)); 610 } 611 x += hx; 612 /*if (max > 7200150000.0) */ 613 /* printf("max very big when i = %d\n",i); */ 614 } 615 PetscCall(PetscDrawAxisSetLabels(axis,"Right hand side","","")); 616 PetscCall(PetscDrawLGSetLegend(lg,NULL)); 617 PetscCall(PetscDrawLGDraw(lg)); 618 619 /* 620 Plot the solution 621 */ 622 PetscCall(PetscDrawLGSetDimension(lg,1)); 623 PetscCall(PetscDrawViewPortsSet(ports,0)); 624 PetscCall(PetscDrawLGReset(lg)); 625 x = hx*xs; 626 PetscCall(PetscDrawLGSetLimits(lg,x,x+(xm-1)*hx,-1.1,1.1)); 627 PetscCall(PetscDrawLGSetColors(lg,colors)); 628 for (i=xs; i<xs+xm; i++) { 629 xx[0] = x; 630 yy[0] = PetscRealPart(u[i]); 631 PetscCall(PetscDrawLGAddPoint(lg,xx,yy)); 632 x += hx; 633 } 634 PetscCall(PetscDrawAxisSetLabels(axis,"Solution","","")); 635 PetscCall(PetscDrawLGDraw(lg)); 636 637 /* 638 Print the forces as arrows on the solution 639 */ 640 x = hx*xs; 641 cnt = xm/60; 642 cnt = (!cnt) ? 1 : cnt; 643 644 for (i=xs; i<xs+xm; i += cnt) { 645 y = yup = ydown = PetscRealPart(u[i]); 646 c = (u[i-1] + u[i+1] - 2.0*u[i])*sx; 647 r = (u[i] + u[i+2] - 2.0*u[i+1])*sx; 648 l = (u[i-2] + u[i] - 2.0*u[i-1])*sx; 649 len = -.5*PetscRealPart(ctx->kappa*(l + r - 2.0*c)*sx)/max; 650 PetscCall(PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_RED)); 651 if (ctx->cahnhillard) { 652 if (len < 0.) ydown += len; 653 else yup += len; 654 655 switch (ctx->energy) { 656 case 1: /* double well */ 657 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; 658 break; 659 case 2: /* double obstacle */ 660 len = -.5*PetscRealPart(u[i-1] + u[i+1] - 2.0*u[i])*sx/max; 661 break; 662 case 3: /* logarithmic + double well */ 663 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; 664 if (len < 0.) ydown += len; 665 else yup += len; 666 667 if (ctx->truncation==2) { /* quadratic */ 668 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; 669 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; 670 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); 671 } else { /* cubic */ 672 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 673 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 674 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); 675 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); 676 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); 677 } 678 y2 = len < 0 ? ydown : yup; 679 PetscCall(PetscDrawArrow(draw,x,y2,x,y2+len2,PETSC_DRAW_PLUM)); 680 break; 681 case 4: /* logarithmic + double obstacle */ 682 len = -.5*theta_c*PetscRealPart(-(u[i-1] + u[i+1] - 2.0*u[i])*sx/max); 683 if (len < 0.) ydown += len; 684 else yup += len; 685 686 if (ctx->truncation==2) { /* quadratic */ 687 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; 688 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; 689 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); 690 } else { /* cubic */ 691 a = 2.0*theta*(1.0-2.0*tol)/(16.0*tol*tol*(1.0-tol)*(1.0-tol)); 692 b = theta/(4.0*tol*(1.0-tol)) - a*(1.0-2.0*tol); 693 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; 694 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; 695 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; 696 } 697 y2 = len < 0 ? ydown : yup; 698 PetscCall(PetscDrawArrow(draw,x,y2,x,y2+len2,PETSC_DRAW_PLUM)); 699 break; 700 } 701 PetscCall(PetscDrawArrow(draw,x,y,x,y+len,PETSC_DRAW_BLUE)); 702 } 703 x += cnt*hx; 704 } 705 PetscCall(DMDAVecRestoreArrayRead(da,localU,&x)); 706 PetscCall(DMRestoreLocalVector(da,&localU)); 707 PetscCall(PetscDrawStringSetSize(draw,.2,.2)); 708 PetscCall(PetscDrawFlush(draw)); 709 PetscCall(PetscDrawSetPause(draw,pause)); 710 PetscCall(PetscDrawPause(draw)); 711 PetscFunctionReturn(0); 712 } 713 714 PetscErrorCode MyDestroy(void **ptr) 715 { 716 UserCtx *ctx = *(UserCtx**)ptr; 717 718 PetscFunctionBegin; 719 PetscCall(PetscDrawViewPortsDestroy(ctx->ports)); 720 PetscFunctionReturn(0); 721 } 722 723 /*TEST 724 725 test: 726 TODO: currently requires initial condition file generated by heat 727 728 TEST*/ 729