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