1 static char help[] = "Solves biharmonic equation in 1d.\n"; 2 3 /* 4 Solves the equation 5 6 u_t = - kappa \Delta \Delta u 7 Periodic boundary conditions 8 9 Evolve the biharmonic heat equation: 10 --------------- 11 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor 12 13 Evolve with the restriction that -1 <= u <= 1; i.e. as a variational inequality 14 --------------- 15 ./biharmonic -ts_monitor -snes_monitor -pc_type lu -draw_pause .1 -snes_converged_reason -draw_pause -2 -ts_type cn -da_refine 5 -mymonitor 16 17 u_t = kappa \Delta \Delta u + 6.*u*(u_x)^2 + (3*u^2 - 12) \Delta u 18 -1 <= u <= 1 19 Periodic boundary conditions 20 21 Evolve the Cahn-Hillard equations: double well Initial hump shrinks then grows 22 --------------- 23 ./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 24 25 Initial hump neither shrinks nor grows when degenerate (otherwise similar solution) 26 27 ./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 28 29 ./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 30 31 Evolve the Cahn-Hillard equations: double obstacle 32 --------------- 33 ./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 34 35 Evolve the Cahn-Hillard equations: logarithmic + double well (never shrinks and then grows) 36 --------------- 37 ./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 38 39 ./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 40 41 Evolve the Cahn-Hillard equations: logarithmic + double obstacle (never shrinks, never grows) 42 --------------- 43 ./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 44 45 */ 46 #include <petscdm.h> 47 #include <petscdmda.h> 48 #include <petscts.h> 49 #include <petscdraw.h> 50 51 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 *); 52 typedef struct { 53 PetscBool cahnhillard; 54 PetscBool degenerate; 55 PetscReal kappa; 56 PetscInt energy; 57 PetscReal tol; 58 PetscReal theta, theta_c; 59 PetscInt truncation; 60 PetscBool netforce; 61 PetscDrawViewPorts *ports; 62 } UserCtx; 63 64 int main(int argc, char **argv) 65 { 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, NULL, 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 { 198 DM da; 199 PetscInt i, Mx, xs, xm; 200 PetscReal hx, sx; 201 PetscScalar *x, *f, c, r, l; 202 Vec localX; 203 UserCtx *ctx = (UserCtx *)ptr; 204 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 */ 205 206 PetscFunctionBegin; 207 PetscCall(TSGetDM(ts, &da)); 208 PetscCall(DMGetLocalVector(da, &localX)); 209 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)); 210 211 hx = 1.0 / (PetscReal)Mx; 212 sx = 1.0 / (hx * hx); 213 214 /* 215 Scatter ghost points to local vector,using the 2-step process 216 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 217 By placing code between these two statements, computations can be 218 done while messages are in transition. 219 */ 220 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 221 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 222 223 /* 224 Get pointers to vector data 225 */ 226 PetscCall(DMDAVecGetArrayRead(da, localX, &x)); 227 PetscCall(DMDAVecGetArray(da, F, &f)); 228 229 /* 230 Get local grid boundaries 231 */ 232 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 233 234 /* 235 Compute function over the locally owned part of the grid 236 */ 237 for (i = xs; i < xs + xm; i++) { 238 if (ctx->degenerate) { 239 c = (1. - x[i] * x[i]) * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx; 240 r = (1. - x[i + 1] * x[i + 1]) * (x[i] + x[i + 2] - 2.0 * x[i + 1]) * sx; 241 l = (1. - x[i - 1] * x[i - 1]) * (x[i - 2] + x[i] - 2.0 * x[i - 1]) * sx; 242 } else { 243 c = (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx; 244 r = (x[i] + x[i + 2] - 2.0 * x[i + 1]) * sx; 245 l = (x[i - 2] + x[i] - 2.0 * x[i - 1]) * sx; 246 } 247 f[i] = -ctx->kappa * (l + r - 2.0 * c) * sx; 248 if (ctx->cahnhillard) { 249 switch (ctx->energy) { 250 case 1: /* double well */ 251 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; 252 break; 253 case 2: /* double obstacle */ 254 f[i] += -(x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx; 255 break; 256 case 3: /* logarithmic + double well */ 257 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; 258 if (ctx->truncation == 2) { /* log function with approximated with a quadratic polynomial outside -1.0+2*tol, 1.0-2*tol */ 259 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; 260 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; 261 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; 262 } else { /* log function is approximated with a cubic polynomial outside -1.0+2*tol, 1.0-2*tol */ 263 a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol)); 264 b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol); 265 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; 266 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; 267 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; 268 } 269 break; 270 case 4: /* logarithmic + double obstacle */ 271 f[i] += -theta_c * (x[i - 1] + x[i + 1] - 2.0 * x[i]) * sx; 272 if (ctx->truncation == 2) { /* quadratic */ 273 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; 274 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; 275 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; 276 } else { /* cubic */ 277 a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol)); 278 b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol); 279 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; 280 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; 281 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; 282 } 283 break; 284 } 285 } 286 } 287 288 /* 289 Restore vectors 290 */ 291 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x)); 292 PetscCall(DMDAVecRestoreArray(da, F, &f)); 293 PetscCall(DMRestoreLocalVector(da, &localX)); 294 PetscFunctionReturn(PETSC_SUCCESS); 295 } 296 297 /* ------------------------------------------------------------------- */ 298 /* 299 FormJacobian - Evaluates nonlinear function's Jacobian 300 301 */ 302 PetscErrorCode FormJacobian(TS ts, PetscReal ftime, Vec X, Mat A, Mat B, void *ptr) 303 { 304 DM da; 305 PetscInt i, Mx, xs, xm; 306 MatStencil row, cols[5]; 307 PetscReal hx, sx; 308 PetscScalar *x, vals[5]; 309 Vec localX; 310 UserCtx *ctx = (UserCtx *)ptr; 311 312 PetscFunctionBegin; 313 PetscCall(TSGetDM(ts, &da)); 314 PetscCall(DMGetLocalVector(da, &localX)); 315 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)); 316 317 hx = 1.0 / (PetscReal)Mx; 318 sx = 1.0 / (hx * hx); 319 320 /* 321 Scatter ghost points to local vector,using the 2-step process 322 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 323 By placing code between these two statements, computations can be 324 done while messages are in transition. 325 */ 326 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 327 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 328 329 /* 330 Get pointers to vector data 331 */ 332 PetscCall(DMDAVecGetArrayRead(da, localX, &x)); 333 334 /* 335 Get local grid boundaries 336 */ 337 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 338 339 /* 340 Compute function over the locally owned part of the grid 341 */ 342 for (i = xs; i < xs + xm; i++) { 343 row.i = i; 344 if (ctx->degenerate) { 345 /*PetscScalar c,r,l; 346 c = (1. - x[i]*x[i])*(x[i-1] + x[i+1] - 2.0*x[i])*sx; 347 r = (1. - x[i+1]*x[i+1])*(x[i] + x[i+2] - 2.0*x[i+1])*sx; 348 l = (1. - x[i-1]*x[i-1])*(x[i-2] + x[i] - 2.0*x[i-1])*sx; */ 349 } else { 350 cols[0].i = i - 2; 351 vals[0] = -ctx->kappa * sx * sx; 352 cols[1].i = i - 1; 353 vals[1] = 4.0 * ctx->kappa * sx * sx; 354 cols[2].i = i; 355 vals[2] = -6.0 * ctx->kappa * sx * sx; 356 cols[3].i = i + 1; 357 vals[3] = 4.0 * ctx->kappa * sx * sx; 358 cols[4].i = i + 2; 359 vals[4] = -ctx->kappa * sx * sx; 360 } 361 PetscCall(MatSetValuesStencil(B, 1, &row, 5, cols, vals, INSERT_VALUES)); 362 363 if (ctx->cahnhillard) { 364 switch (ctx->energy) { 365 case 1: /* double well */ 366 /* 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; */ 367 break; 368 case 2: /* double obstacle */ 369 /* f[i] += -(x[i-1] + x[i+1] - 2.0*x[i])*sx; */ 370 break; 371 case 3: /* logarithmic + double well */ 372 break; 373 case 4: /* logarithmic + double obstacle */ 374 break; 375 } 376 } 377 } 378 379 /* 380 Restore vectors 381 */ 382 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x)); 383 PetscCall(DMRestoreLocalVector(da, &localX)); 384 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 385 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 386 if (A != B) { 387 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 388 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 389 } 390 PetscFunctionReturn(PETSC_SUCCESS); 391 } 392 /* ------------------------------------------------------------------- */ 393 PetscErrorCode FormInitialSolution(DM da, Vec U) 394 { 395 PetscInt i, xs, xm, Mx, N, scale; 396 PetscScalar *u; 397 PetscReal r, hx, x; 398 const PetscScalar *f; 399 Vec finesolution; 400 PetscViewer viewer; 401 402 PetscFunctionBegin; 403 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)); 404 405 hx = 1.0 / (PetscReal)Mx; 406 407 /* 408 Get pointers to vector data 409 */ 410 PetscCall(DMDAVecGetArray(da, U, &u)); 411 412 /* 413 Get local grid boundaries 414 */ 415 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 416 417 /* 418 Seee heat.c for how to generate InitialSolution.heat 419 */ 420 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "InitialSolution.heat", FILE_MODE_READ, &viewer)); 421 PetscCall(VecCreate(PETSC_COMM_WORLD, &finesolution)); 422 PetscCall(VecLoad(finesolution, viewer)); 423 PetscCall(PetscViewerDestroy(&viewer)); 424 PetscCall(VecGetSize(finesolution, &N)); 425 scale = N / Mx; 426 PetscCall(VecGetArrayRead(finesolution, &f)); 427 428 /* 429 Compute function over the locally owned part of the grid 430 */ 431 for (i = xs; i < xs + xm; i++) { 432 x = i * hx; 433 r = PetscSqrtReal((x - .5) * (x - .5)); 434 if (r < .125) u[i] = 1.0; 435 else u[i] = -.5; 436 437 /* With the initial condition above the method is first order in space */ 438 /* this is a smooth initial condition so the method becomes second order in space */ 439 /*u[i] = PetscSinScalar(2*PETSC_PI*x); */ 440 u[i] = f[scale * i]; 441 } 442 PetscCall(VecRestoreArrayRead(finesolution, &f)); 443 PetscCall(VecDestroy(&finesolution)); 444 445 /* 446 Restore vectors 447 */ 448 PetscCall(DMDAVecRestoreArray(da, U, &u)); 449 PetscFunctionReturn(PETSC_SUCCESS); 450 } 451 452 /* 453 This routine is not parallel 454 */ 455 PetscErrorCode MyMonitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ptr) 456 { 457 UserCtx *ctx = (UserCtx *)ptr; 458 PetscDrawLG lg; 459 PetscScalar *u, l, r, c; 460 PetscInt Mx, i, xs, xm, cnt; 461 PetscReal x, y, hx, pause, sx, len, max, xx[4], yy[4], xx_netforce, yy_netforce, yup, ydown, y2, len2; 462 PetscDraw draw; 463 Vec localU; 464 DM da; 465 int colors[] = {PETSC_DRAW_YELLOW, PETSC_DRAW_RED, PETSC_DRAW_BLUE, PETSC_DRAW_PLUM, PETSC_DRAW_BLACK}; 466 /* 467 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"}}; 468 */ 469 PetscDrawAxis axis; 470 PetscDrawViewPorts *ports; 471 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 */ 472 PetscReal vbounds[] = {-1.1, 1.1}; 473 474 PetscFunctionBegin; 475 PetscCall(PetscViewerDrawSetBounds(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, vbounds)); 476 PetscCall(PetscViewerDrawResize(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 800, 600)); 477 PetscCall(TSGetDM(ts, &da)); 478 PetscCall(DMGetLocalVector(da, &localU)); 479 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)); 480 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL)); 481 hx = 1.0 / (PetscReal)Mx; 482 sx = 1.0 / (hx * hx); 483 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU)); 484 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU)); 485 PetscCall(DMDAVecGetArrayRead(da, localU, &u)); 486 487 PetscCall(PetscViewerDrawGetDrawLG(PETSC_VIEWER_DRAW_(PETSC_COMM_WORLD), 1, &lg)); 488 PetscCall(PetscDrawLGGetDraw(lg, &draw)); 489 PetscCall(PetscDrawCheckResizedWindow(draw)); 490 if (!ctx->ports) PetscCall(PetscDrawViewPortsCreateRect(draw, 1, 3, &ctx->ports)); 491 ports = ctx->ports; 492 PetscCall(PetscDrawLGGetAxis(lg, &axis)); 493 PetscCall(PetscDrawLGReset(lg)); 494 495 xx[0] = 0.0; 496 xx[1] = 1.0; 497 cnt = 2; 498 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-zoom", xx, &cnt, NULL)); 499 xs = xx[0] / hx; 500 xm = (xx[1] - xx[0]) / hx; 501 502 /* 503 Plot the energies 504 */ 505 PetscCall(PetscDrawLGSetDimension(lg, 1 + (ctx->cahnhillard ? 1 : 0) + (ctx->energy == 3))); 506 PetscCall(PetscDrawLGSetColors(lg, colors + 1)); 507 PetscCall(PetscDrawViewPortsSet(ports, 2)); 508 x = hx * xs; 509 for (i = xs; i < xs + xm; i++) { 510 xx[0] = xx[1] = xx[2] = x; 511 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); 512 else yy[0] = PetscRealPart(.25 * ctx->kappa * (u[i - 1] - u[i + 1]) * (u[i - 1] - u[i + 1]) * sx); 513 514 if (ctx->cahnhillard) { 515 switch (ctx->energy) { 516 case 1: /* double well */ 517 yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - u[i] * u[i])); 518 break; 519 case 2: /* double obstacle */ 520 yy[1] = .5 * PetscRealPart(1. - u[i] * u[i]); 521 break; 522 case 3: /* logarithm + double well */ 523 yy[1] = .25 * PetscRealPart((1. - u[i] * u[i]) * (1. - 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 case 4: /* logarithm + double obstacle */ 529 yy[1] = .5 * theta_c * PetscRealPart(1.0 - u[i] * u[i]); 530 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)); 531 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)); 532 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)); 533 break; 534 default: 535 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "It will always be one of the values"); 536 } 537 } 538 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 539 x += hx; 540 } 541 PetscCall(PetscDrawGetPause(draw, &pause)); 542 PetscCall(PetscDrawSetPause(draw, 0.0)); 543 PetscCall(PetscDrawAxisSetLabels(axis, "Energy", "", "")); 544 /* PetscCall(PetscDrawLGSetLegend(lg,legend[ctx->energy-1])); */ 545 PetscCall(PetscDrawLGDraw(lg)); 546 547 /* 548 Plot the forces 549 */ 550 PetscCall(PetscDrawLGSetDimension(lg, 0 + (ctx->cahnhillard ? 2 : 0) + (ctx->energy == 3))); 551 PetscCall(PetscDrawLGSetColors(lg, colors + 1)); 552 PetscCall(PetscDrawViewPortsSet(ports, 1)); 553 PetscCall(PetscDrawLGReset(lg)); 554 x = xs * hx; 555 max = 0.; 556 for (i = xs; i < xs + xm; i++) { 557 xx[0] = xx[1] = xx[2] = xx[3] = x; 558 xx_netforce = x; 559 if (ctx->degenerate) { 560 c = (1. - u[i] * u[i]) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx; 561 r = (1. - u[i + 1] * u[i + 1]) * (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx; 562 l = (1. - u[i - 1] * u[i - 1]) * (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx; 563 } else { 564 c = (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx; 565 r = (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx; 566 l = (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx; 567 } 568 yy[0] = PetscRealPart(-ctx->kappa * (l + r - 2.0 * c) * sx); 569 yy_netforce = yy[0]; 570 max = PetscMax(max, PetscAbs(yy[0])); 571 if (ctx->cahnhillard) { 572 switch (ctx->energy) { 573 case 1: /* double well */ 574 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); 575 break; 576 case 2: /* double obstacle */ 577 yy[1] = -PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx; 578 break; 579 case 3: /* logarithmic + double well */ 580 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); 581 if (ctx->truncation == 2) { /* quadratic */ 582 if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx; 583 else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = (.25 * theta / (tol - tol * tol)) * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx; 584 else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx); 585 } else { /* cubic */ 586 a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol)); 587 b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol); 588 if (PetscRealPart(u[i]) < -1.0 + 2.0 * tol) yy[2] = PetscRealPart(-1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (-1.0 * a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx); 589 else if (PetscRealPart(u[i]) > 1.0 - 2.0 * tol) yy[2] = PetscRealPart(1.0 * a * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (a * u[i] + b) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx); 590 else yy[2] = PetscRealPart(2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx); 591 } 592 break; 593 case 4: /* logarithmic + double obstacle */ 594 yy[1] = theta_c * PetscRealPart(-(u[i - 1] + u[i + 1] - 2.0 * u[i])) * sx; 595 if (ctx->truncation == 2) { 596 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; 597 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; 598 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); 599 } else { 600 a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol)); 601 b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol); 602 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); 603 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); 604 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); 605 } 606 break; 607 default: 608 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "It will always be one of the values"); 609 } 610 if (ctx->energy < 3) { 611 max = PetscMax(max, PetscAbs(yy[1])); 612 yy[2] = yy[0] + yy[1]; 613 yy_netforce = yy[2]; 614 } else { 615 max = PetscMax(max, PetscAbs(yy[1] + yy[2])); 616 yy[3] = yy[0] + yy[1] + yy[2]; 617 yy_netforce = yy[3]; 618 } 619 } 620 if (ctx->netforce) { 621 PetscCall(PetscDrawLGAddPoint(lg, &xx_netforce, &yy_netforce)); 622 } else { 623 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 624 } 625 x += hx; 626 /*if (max > 7200150000.0) */ 627 /* printf("max very big when i = %d\n",i); */ 628 } 629 PetscCall(PetscDrawAxisSetLabels(axis, "Right hand side", "", "")); 630 PetscCall(PetscDrawLGSetLegend(lg, NULL)); 631 PetscCall(PetscDrawLGDraw(lg)); 632 633 /* 634 Plot the solution 635 */ 636 PetscCall(PetscDrawLGSetDimension(lg, 1)); 637 PetscCall(PetscDrawViewPortsSet(ports, 0)); 638 PetscCall(PetscDrawLGReset(lg)); 639 x = hx * xs; 640 PetscCall(PetscDrawLGSetLimits(lg, x, x + (xm - 1) * hx, -1.1, 1.1)); 641 PetscCall(PetscDrawLGSetColors(lg, colors)); 642 for (i = xs; i < xs + xm; i++) { 643 xx[0] = x; 644 yy[0] = PetscRealPart(u[i]); 645 PetscCall(PetscDrawLGAddPoint(lg, xx, yy)); 646 x += hx; 647 } 648 PetscCall(PetscDrawAxisSetLabels(axis, "Solution", "", "")); 649 PetscCall(PetscDrawLGDraw(lg)); 650 651 /* 652 Print the forces as arrows on the solution 653 */ 654 x = hx * xs; 655 cnt = xm / 60; 656 cnt = (!cnt) ? 1 : cnt; 657 658 for (i = xs; i < xs + xm; i += cnt) { 659 y = yup = ydown = PetscRealPart(u[i]); 660 c = (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx; 661 r = (u[i] + u[i + 2] - 2.0 * u[i + 1]) * sx; 662 l = (u[i - 2] + u[i] - 2.0 * u[i - 1]) * sx; 663 len = -.5 * PetscRealPart(ctx->kappa * (l + r - 2.0 * c) * sx) / max; 664 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_RED)); 665 if (ctx->cahnhillard) { 666 if (len < 0.) ydown += len; 667 else yup += len; 668 669 switch (ctx->energy) { 670 case 1: /* double well */ 671 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; 672 break; 673 case 2: /* double obstacle */ 674 len = -.5 * PetscRealPart(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max; 675 break; 676 case 3: /* logarithmic + double well */ 677 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; 678 if (len < 0.) ydown += len; 679 else yup += len; 680 681 if (ctx->truncation == 2) { /* quadratic */ 682 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; 683 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; 684 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); 685 } else { /* cubic */ 686 a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol)); 687 b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol); 688 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); 689 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); 690 else len2 = PetscRealPart(.5 * (2.0 * theta * u[i] / ((1.0 - u[i] * u[i]) * (1.0 - u[i] * u[i])) * .25 * (u[i + 1] - u[i - 1]) * (u[i + 1] - u[i - 1]) * sx + (theta / (1.0 - u[i] * u[i])) * (u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx) / max); 691 } 692 y2 = len < 0 ? ydown : yup; 693 PetscCall(PetscDrawArrow(draw, x, y2, x, y2 + len2, PETSC_DRAW_PLUM)); 694 break; 695 case 4: /* logarithmic + double obstacle */ 696 len = -.5 * theta_c * PetscRealPart(-(u[i - 1] + u[i + 1] - 2.0 * u[i]) * sx / max); 697 if (len < 0.) ydown += len; 698 else yup += len; 699 700 if (ctx->truncation == 2) { /* quadratic */ 701 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; 702 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; 703 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); 704 } else { /* cubic */ 705 a = 2.0 * theta * (1.0 - 2.0 * tol) / (16.0 * tol * tol * (1.0 - tol) * (1.0 - tol)); 706 b = theta / (4.0 * tol * (1.0 - tol)) - a * (1.0 - 2.0 * tol); 707 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; 708 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; 709 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; 710 } 711 y2 = len < 0 ? ydown : yup; 712 PetscCall(PetscDrawArrow(draw, x, y2, x, y2 + len2, PETSC_DRAW_PLUM)); 713 break; 714 } 715 PetscCall(PetscDrawArrow(draw, x, y, x, y + len, PETSC_DRAW_BLUE)); 716 } 717 x += cnt * hx; 718 } 719 PetscCall(DMDAVecRestoreArrayRead(da, localU, &x)); 720 PetscCall(DMRestoreLocalVector(da, &localU)); 721 PetscCall(PetscDrawStringSetSize(draw, .2, .2)); 722 PetscCall(PetscDrawFlush(draw)); 723 PetscCall(PetscDrawSetPause(draw, pause)); 724 PetscCall(PetscDrawPause(draw)); 725 PetscFunctionReturn(PETSC_SUCCESS); 726 } 727 728 PetscErrorCode MyDestroy(void **ptr) 729 { 730 UserCtx *ctx = *(UserCtx **)ptr; 731 732 PetscFunctionBegin; 733 PetscCall(PetscDrawViewPortsDestroy(ctx->ports)); 734 PetscFunctionReturn(PETSC_SUCCESS); 735 } 736 737 /*TEST 738 739 test: 740 TODO: currently requires initial condition file generated by heat 741 742 TEST*/ 743