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