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