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