1 | /*************************************** 2 | $Header$ 3 | 4 | This file contains the heart of the Cahn-Hilliard formulation, in particular 5 | the functions which build the finite difference residuals and Jacobian. 6 | ***************************************/ 7 | 8 | 9 | /* Including cahnhill.h includes the necessary PETSc header files */ 10 | #include "cahnhill.h" 11 | 12 | 13 | /*++++++++++++++++++++++++++++++++++++++ 14 | This calculates the derivative of homogeneous free energy with respect to 15 | concentration, which is a component of the chemical potential. It currently 16 | uses 17 | +latex+$$\Psi' = C(1-C)\left(\frac{1}{2}+m-C\right)$$ 18 | +html+ <center>Psi' = <i>C</i> (1-<i>C</i>) (0.5+<i>m</i>+<i>C</i>)</center> 19 | which gives (meta)stable equilibria at 20 | +latex+$C=0$ and 1 and an unstable equilibrium at $C=\frac{1}{2}+m$; if $m>0$ 21 | +html+ <i>C</i>=0 and 1 and an unstable equilibrium at <i>C</i> = 1/2 + 22 | +html+ <i>m</i>; if <i>m</i>>0 23 | then the 0 phase is stable and vice versa. 24 | 25 | PetscScalar ch_psiprime It returns the derivative itself. 26 | 27 | PetscScalar C The concentration. 28 | 29 | PetscScalar mparam The model parameter 30 | +latex+$m$. 31 | +html+ <i>m</i>. 32 | ++++++++++++++++++++++++++++++++++++++*/ 33 | 34 | static inline PetscScalar ch_psiprime (PetscScalar C, PetscScalar mparam) 35 | { return C*(1-C)*(0.5+mparam-C); } 36 | 37 | #define PSIPRIME_FLOPS 5 38 | 39 | 40 | /*++++++++++++++++++++++++++++++++++++++ 41 | This calculates the second derivative of homogeneous free energy with respect 42 | to concentration, for insertion into the Jacobian matrix. See the 43 | +latex+{\tt ch\_psiprime()} function for the definition of $\Psi$. 44 | +html+ <tt>ch_psiprime()</tt> function for the definition of Psi. 45 | 46 | PetscScalar ch_psidoubleprime It returns the second derivative 47 | +latex+$\Psi''(C)$. 48 | -latex-Psi''(C). 49 | 50 | PetscScalar C The concentration. 51 | 52 | PetscScalar mparam The model parameter 53 | +latex+$m$. 54 | +html+ <i>m</i>. 55 | ++++++++++++++++++++++++++++++++++++++*/ 56 | 57 | static inline PetscScalar ch_psidoubleprime (PetscScalar C, PetscScalar mparam) 58 | { return 3*C*C - (3+2*mparam)*C + 0.5+mparam; } 59 | 60 | #define PSIDOUBLEPRIME_FLOPS 8 61 | 62 | 63 | /*++++++++++++++++++++++++++++++++++++++ 64 | This function computes the residual from indices to points in the 65 | concentration array. ``Up'' refers to the positive 66 | +latex+$y$-direction, ``down'' to negative $y$, ``left'' to negative $x$ and 67 | +latex+``right'' to positive $x$. 68 | +html+ <i>y</i>-direction, ``down'' to negative <i>y</i>, ``left'' to 69 | +html+ negative <i>x</i> and ``right'' to positive <i>x</i>. 70 | 71 | PetscScalar ch_residual_2d Returns the residual itself 72 | 73 | PetscScalar *conc Array of concentrations 74 | 75 | PetscScalar alpha Model parameter 76 | +latex+$\kappa\alpha$ 77 | -latex-kappa alpha 78 | 79 | PetscScalar beta Model parameter 80 | +latex+$\kappa\beta$ 81 | -latex-kappa beta 82 | 83 | PetscScalar mparam Model parameter 84 | +latex+$m$ 85 | -latex-m 86 | 87 | PetscScalar hx Inverse square 88 | +latex+$x$-spacing $h_x^{-2}$ 89 | +html+ <i>x</i>-spacing <i>h<sub>x</sub></i><sup>-2</sup> 90 | 91 | PetscScalar hy Inverse square 92 | +latex+$y$-spacing $h_y^{-2}$ 93 | +html+ <i>y</i>-spacing <i>h<sub>y</sub></i><sup>-2</sup> 94 | 95 | int upup Index to array position two cells up from current 96 | 97 | int upleft Index to array position one cell up and one left from current 98 | 99 | int up Index to array position one cell up from current 100 | 101 | int upright Index to array position one cell up and one right from current 102 | 103 | int leftleft Index to array position two cells left of current 104 | 105 | int left Index to array position one cell left of current 106 | 107 | int current Index to current cell array position 108 | 109 | int right Index to array position one cell right of current 110 | 111 | int rightright Index to array position two cells right of current 112 | 113 | int downleft Index to array position one cell down and one left from current 114 | 115 | int down Index to array position one cell down from current 116 | 117 | int downright Index to array position one cell down and one right from current 118 | 119 | int downdown Index to array position two cells down from current 120 | ++++++++++++++++++++++++++++++++++++++*/ 121 | 122 | static inline PetscScalar ch_residual_2d 123 | (PetscScalar *conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam, 124 | PetscScalar hx, PetscScalar hy, int upup, int upleft, int up, int upright, int 125 | leftleft, int left, int current, int right, int rightright, int downleft, int 126 | down, int downright, int downdown) 127 | { 128 | PetscScalar retval, hx2, hy2, hxhy; 129 | 130 | hx2 = hx*hx; hy2 = hy*hy; hxhy = hx*hy; 131 | 132 | /*+ This calculates the 133 | +latex+$\beta$-term, $\kappa\beta$ times the Laplacian of $\Psi'(C)$, 134 | -latex-beta term, kappa*beta times the laplacian of Psi prime(C), 135 | +*/ 136 | retval = beta * 137 | (hx * (ch_psiprime(conc[left], mparam) + ch_psiprime(conc[right], mparam)) 138 | + hy * (ch_psiprime(conc[up], mparam) + ch_psiprime(conc[down], mparam)) 139 | - 2*(hx+hy) * ch_psiprime(conc[current], mparam)); 140 | 141 | /*+ then subtracts the 142 | +latex+$\alpha$-term, $\kappa\alpha\nabla^2\nabla^2C$. 143 | -latex-alpha term, kappa*alpha times the Laplacian of the Laplacian of C. 144 | +*/ 145 | retval -= alpha * 146 | (hx2 * (conc[leftleft] + conc[rightright]) 147 | + hy2 * (conc[upup] + conc[downdown]) 148 | + 2*hxhy * (conc[upleft]+conc[upright]+conc[downleft]+conc[downright]) 149 | - (4*hx2 + 4*hxhy) * (conc[left] + conc[right]) 150 | - (4*hy2 + 4*hxhy) * (conc[up] + conc[down]) 151 | + (6*hx2 + 6*hy2 + 8*hxhy) * conc[current]); 152 | 153 | return retval; 154 | } 155 | 156 | #define RESIDUAL_FLOPS_2D (5*PSIPRIME_FLOPS + 10 + 32) 157 | 158 | 159 | /*++++++++++++++++++++++++++++++++++++++ 160 | This function computes the residual from indices to points in the 161 | concentration array. ``Front refers to the positive 162 | +latex+$z$-direction, ``back'' to negative $z$, ``up'' to positive 163 | +latex+$y$, ``down'' to negative $y$, ``left'' to 164 | +latex+negative $x$ and ``right'' to positive $x$. 165 | +html+ <i>z</i>-direction, ``back'' to negative <i>z</i>, ``up'' to positive 166 | +html+ <i>y</i>, ``down'' to negative <i>y</i>, ``left'' to 167 | +html+ negative <i>x</i> and ``right'' to positive <i>x</i>. 168 | 169 | PetscScalar ch_residual_3d Returns the residual itself 170 | 171 | PetscScalar *conc Array of concentrations 172 | 173 | PetscScalar alpha Model parameter 174 | +latex+$\kappa\alpha$ 175 | -latex-kappa alpha 176 | 177 | PetscScalar beta Model parameter 178 | +latex+$\kappa\beta$ 179 | -latex-kappa beta 180 | 181 | PetscScalar mparam Model parameter 182 | +latex+$m$ 183 | -latex-m 184 | 185 | PetscScalar hx Inverse square 186 | +latex+$x$-spacing $h_x^{-2}$ 187 | +html+ <i>x</i>-spacing <i>h<sub>x</sub></i><sup>-2</sup> 188 | 189 | PetscScalar hy Inverse square 190 | +latex+$y$-spacing $h_y^{-2}$ 191 | +html+ <i>y</i>-spacing <i>h<sub>y</sub></i><sup>-2</sup> 192 | 193 | PetscScalar hz Inverse square 194 | +latex+$z$-spacing $h_z^{-2}$ 195 | +html+ <i>z</i>-spacing <i>h<sub>z</sub></i><sup>-2</sup> 196 | 197 | int frontfront Index to array position two cells in front of current 198 | 199 | int frontup Index to array position one cell front and one up from current 200 | 201 | int frontleft Index to array position one cell front and one left from current 202 | 203 | int front Index to array position one cell in front of current 204 | 205 | int frontright Index to array position one cell front and one right from current 206 | 207 | int frontdown Index to array position one cell front and one down from current 208 | 209 | int upup Index to array position two cells up from current 210 | 211 | int upleft Index to array position one cell up and one left from current 212 | 213 | int up Index to array position one cell up from current 214 | 215 | int upright Index to array position one cell up and one right from current 216 | 217 | int leftleft Index to array position two cells left of current 218 | 219 | int left Index to array position one cell left of current 220 | 221 | int current Index to current cell array position 222 | 223 | int right Index to array position one cell right of current 224 | 225 | int rightright Index to array position two cells right of current 226 | 227 | int downleft Index to array position one cell down and one left from current 228 | 229 | int down Index to array position one cell down from current 230 | 231 | int downright Index to array position one cell down and one right from current 232 | 233 | int downdown Index to array position two cells down from current 234 | 235 | int backup Index to array position one cell back and one up from current 236 | 237 | int backleft Index to array position one cell back and one left from current 238 | 239 | int back Index to array position one cell in back of current 240 | 241 | int backright Index to array position one cell back and one right from current 242 | 243 | int backdown Index to array position one cell back and one down from current 244 | 245 | int backback Index to array position two cells in back of current 246 | ++++++++++++++++++++++++++++++++++++++*/ 247 | 248 | static inline PetscScalar ch_residual_3d 249 | (PetscScalar *conc, PetscScalar alpha, PetscScalar beta, PetscScalar mparam, 250 | PetscScalar hx, PetscScalar hy, PetscScalar hz, int frontfront, 251 | int frontup, int frontleft, int front, int frontright, int frontdown, 252 | int upup, int upleft, int up, int upright, int leftleft, int left, 253 | int current, int right, int rightright, int downleft, int down, int downright, 254 | int downdown, 255 | int backup, int backleft, int back, int backright, int backdown, 256 | int backback) 257 | { 258 | PetscScalar retval, hx2, hy2, hz2, hxhy, hxhz, hyhz; 259 | 260 | hx2 = hx*hx; hy2 = hy*hy; hz2 = hz*hz; 261 | hxhy = hx*hy; hxhz = hx*hz; hyhz = hy*hz; 262 | 263 | /*+ This calculates the 264 | +latex+$\beta$-term, $\kappa\beta$ times the Laplacian of $\Psi'(C)$, 265 | -latex-beta term, kappa*beta times the laplacian of Psi prime(C), 266 | +*/ 267 | retval = beta * 268 | (hx * (ch_psiprime(conc[left], mparam) + ch_psiprime(conc[right], mparam)) 269 | + hy * (ch_psiprime(conc[up], mparam) + ch_psiprime(conc[down], mparam)) 270 | + hz * (ch_psiprime(conc[front], mparam)+ ch_psiprime(conc[back], mparam)) 271 | - 2*(hx+hy+hz) * ch_psiprime(conc[current], mparam)); 272 | 273 | /*+ then subtracts the 274 | +latex+$\alpha$-term, $\kappa\alpha\nabla^2\nabla^2C$. 275 | -latex-alpha term, kappa*alpha times the Laplacian of the Laplacian of C. 276 | +*/ 277 | retval -= alpha * 278 | (hx2 * (conc[leftleft] + conc[rightright]) 279 | + hy2 * (conc[upup] + conc[downdown]) 280 | + hz2 * (conc[frontfront] + conc[backback]) 281 | + 2*hxhy * (conc[upleft]+conc[upright]+conc[downleft]+conc[downright]) 282 | + 2*hxhz*(conc[frontleft]+conc[frontright]+conc[backleft]+conc[backright]) 283 | + 2*hyhz * (conc[frontup]+conc[frontdown]+conc[backup]+conc[backdown]) 284 | - (4*hx2 + 4*hxhy + 4*hxhz) * (conc[left] + conc[right]) 285 | - (4*hy2 + 4*hxhy + 4*hyhz) * (conc[up] + conc[down]) 286 | - (4*hz2 + 4*hxhz + 4*hyhz) * (conc[front] + conc[back]) 287 | + (6*hx2 + 6*hy2 + 6*hz2 + 8*hxhy + 8*hxhz + 8*hyhz) * conc[current]); 288 | 289 | return retval; 290 | } 291 | 292 | #define RESIDUAL_FLOPS_3D (7*PSIPRIME_FLOPS + 14 + 65) 293 | 294 | 295 | /*++++++++++++++++++++++++++++++++++++++ 296 | This evaluates the nonlinear finite difference approximation to the residuals 297 | +latex+$R_i$. 298 | +html+ <i>R<sub>i</sub></i>. 299 | Note that it loops on the interior points and the boundary separately, to 300 | avoid conditional statements within the double loop over the local grid 301 | indices. 302 | 303 | int ch_residual_vector_2d It returns zero or an error value. 304 | 305 | Vec X The pre-allocated local vector of unknowns. 306 | 307 | Vec F The pre-allocated local vector of residuals, filled by this function. 308 | 309 | void *ptr Data passed in the application context. 310 | ++++++++++++++++++++++++++++++++++++++*/ 311 | 312 | int ch_residual_vector_2d (Vec X, Vec F, void *ptr) 313 | { 314 | AppCtx *user = (AppCtx*)ptr; 315 | int ierr,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym; 316 | int xints,xinte,yints,yinte; 317 | PetscScalar hx,hy,dhx,dhy,hxm2,hym2; 318 | PetscScalar alpha, beta, mparam; 319 | PetscScalar *x,*f; 320 | Vec localX = user->localX,localF = user->localF; 321 | 322 | /* Scatter ghost points to local vector, using the 2-step process 323 | DAGlobalToLocalBegin(), DAGlobalToLocalEnd(). 324 | By placing code between these two statements, computations can be 325 | done while messages are in transition. */ 326 | ierr = DAGlobalToLocalBegin (user->da, X, INSERT_VALUES, localX); 327 | CHKERRQ (ierr); 328 | 329 | /* Define mesh intervals ratios for uniform grid. 330 | [Note: FD formulae below may someday be normalized by multiplying through 331 | by local volume element to obtain coefficients O(1) in two dimensions.] */ 332 | 333 | mc = user->mc; 334 | chvar = user->chvar; 335 | mx = user->mx; 336 | my = user->my; 337 | mparam = user->mparam; 338 | alpha = user->kappa * user->gamma * user->epsilon; 339 | beta = user->kappa * user->gamma / user->epsilon; 340 | dhx = (PetscScalar)(mx-1); 341 | dhy = (PetscScalar)(my-1); 342 | /* These next two lines hard-code the 1x1 square */ 343 | hx = 1./dhx; 344 | hy = 1./dhy; 345 | hxm2 = 1./hx/hx; 346 | hym2 = 1./hy/hy; 347 | 348 | ierr = DAGlobalToLocalEnd (user->da, X, INSERT_VALUES, localX); 349 | CHKERRQ (ierr); 350 | 351 | /* Get pointers to vector data */ 352 | ierr = VecGetArray (localX, &x); CHKERRQ (ierr); 353 | ierr = VecGetArray (localF, &f); CHKERRQ (ierr); 354 | 355 | /* Get local grid boundaries */ 356 | ierr = DAGetCorners (user->da, &xs, &ys, PETSC_NULL, &xm, &ym, PETSC_NULL); 357 | CHKERRQ (ierr); 358 | ierr = DAGetGhostCorners (user->da, &gxs, &gys, PETSC_NULL, &gxm, &gym, 359 | PETSC_NULL); CHKERRQ (ierr); 360 | 361 | /* Determine the interior of the locally owned part of the grid. */ 362 | if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) { 363 | xints = xs; xinte = xs+xm; } 364 | else { 365 | xints = (xs==0) ? 2:xs; xinte = (xs+xm==mx) ? xs+xm-2 : xs+xm; } 366 | if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) { 367 | yints = ys; yinte = ys+ym; } 368 | else { 369 | yints = (ys==0) ? 2:ys; yinte = (ys+ym==my) ? ys+ym-2 : ys+ym; } 370 | 371 | /* Loop over the interior points */ 372 | for (j=yints-gys; j<yinte-gys; j++) 373 | for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) 374 | f[C(i)] = ch_residual_2d 375 | (x, alpha, beta, mparam, hxm2, hym2, 376 | C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1), 377 | C(i-2), C(i-1), C(i), C(i+1), C(i+2), 378 | C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm)); 379 | 380 | /* Okay, that was the easy part. Now for the fun part. 381 | The following conditionals implement symmetry boundary conditions. */ 382 | 383 | /* Test whether we are on the bottom edge of the global array */ 384 | if (yints == 2) /* Not ys==0 because that would be true for periodic too */ 385 | { 386 | printf ("This must only be reached if we are NOT y-periodic\n"); 387 | /* Bottom two edge lines */ 388 | for (i=xints-gxs; i<xinte-gxs; i++) 389 | { 390 | f[C(i)] = ch_residual_2d 391 | (x, alpha, beta, mparam, hxm2, hym2, 392 | C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1), 393 | C(i-2), C(i-1), C(i), C(i+1), C(i+2), 394 | C(i+gxm-1), C(i+gxm), C(i+gxm+1), C(i+2*gxm)); 395 | f[C(i+gxm)] = ch_residual_2d 396 | (x, alpha, beta, mparam, hxm2, hym2, 397 | C(i+3*gxm), C(i+2*gxm-1), C(i+2*gxm), C(i+2*gxm+1), 398 | C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm+1), C(i+gxm+2), 399 | C(i-1), C(i), C(i+1), C(i+gxm)); 400 | } 401 | 402 | /* Bottom left corner */ 403 | if (xints == 2) 404 | { 405 | f[C(0)] = ch_residual_2d 406 | (x, alpha, beta, mparam, hxm2, hym2, 407 | C(2*gxm), C(gxm+1), C(gxm), C(gxm+1), 408 | C(2), C(1), C(0), C(1), C(2), 409 | C(gxm+1), C(gxm), C(gxm+1), C(2*gxm)); 410 | f[C(1)] = ch_residual_2d 411 | (x, alpha, beta, mparam, hxm2, hym2, 412 | C(2*gxm+1), C(gxm), C(gxm+1), C(gxm+2), 413 | C(1), C(0), C(1), C(2), C(3), 414 | C(gxm), C(gxm+1), C(gxm+2), C(2*gxm+1)); 415 | f[C(gxm)] = ch_residual_2d 416 | (x, alpha, beta, mparam, hxm2, hym2, 417 | C(3*gxm), C(2*gxm+1), C(2*gxm), C(2*gxm+1), 418 | C(gxm+2), C(gxm+1), C(gxm+0), C(gxm+1), C(gxm+2), 419 | C(1), C(0), C(1), C(gxm)); 420 | f[C(gxm+1)] = ch_residual_2d 421 | (x, alpha, beta, mparam, hxm2, hym2, 422 | C(3*gxm+1), C(2*gxm), C(2*gxm+1), C(2*gxm+2), 423 | C(gxm+1), C(gxm), C(gxm+1), C(gxm+2), C(gxm+3), 424 | C(0), C(1), C(2), C(gxm+1)); 425 | } 426 | 427 | /* Bottom right corner */ 428 | if (xinte == mx-2) 429 | { 430 | i = mx-gxs-1; /* Array index of the bottom right corner point */ 431 | f[C(i)] = ch_residual_2d 432 | (x, alpha, beta, mparam, hxm2, hym2, 433 | C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm-1), 434 | C(i-2), C(i-1), C(i), C(i-1), C(i-2), 435 | C(i+gxm-1), C(i+gxm), C(i+gxm-1), C(i+2*gxm)); 436 | f[C(i-1)] = ch_residual_2d 437 | (x, alpha, beta, mparam, hxm2, hym2, 438 | C(i+2*gxm-1), C(i+gxm-2), C(i+gxm-1), C(i+gxm), 439 | C(i-3), C(i-2), C(i-1), C(i), C(i-1), 440 | C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+2*gxm-1)); 441 | f[C(i+gxm)] = ch_residual_2d 442 | (x, alpha, beta, mparam, hxm2, hym2, 443 | C(i+3*gxm), C(i+2*gxm-1), C(i+2*gxm), C(i+2*gxm-1), 444 | C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm-1), C(i+gxm-2), 445 | C(i-1), C(i), C(i-1), C(i+gxm)); 446 | f[C(i+gxm-1)] = ch_residual_2d 447 | (x, alpha, beta, mparam, hxm2, hym2, 448 | C(i+3*gxm-1), C(i+2*gxm-2), C(i+2*gxm-1), C(i+2*gxm), 449 | C(i+gxm-3), C(i+gxm-2), C(i+gxm-1), C(i+gxm), C(i+gxm-1), 450 | C(i-2), C(i-1), C(i), C(i+gxm-1)); 451 | } 452 | } 453 | 454 | /* Test whether we are on the top edge of the global array */ 455 | if (yinte == my-2) 456 | { 457 | printf ("This must only be reached if we are NOT y-periodic\n"); 458 | /* Top two edge lines */ 459 | for (i=(my-gys-1)*gxm + xints-gxs; i<(my-gys-1)*gxm + xinte-gxs; i++) 460 | { 461 | f[C(i)] = ch_residual_2d 462 | (x, alpha, beta, mparam, hxm2, hym2, 463 | C(i-2*gxm), C(i-gxm-1), C(i-gxm), C(i-gxm+1), 464 | C(i-2), C(i-1), C(i), C(i+1), C(i+2), 465 | C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm)); 466 | f[C(i-gxm)] = ch_residual_2d 467 | (x, alpha, beta, mparam, hxm2, hym2, 468 | C(i-gxm), C(i-1), C(i), C(i+1), 469 | C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-gxm+2), 470 | C(i-2*gxm-1), C(i-2*gxm), C(i-2*gxm+1), C(i-3*gxm)); 471 | } 472 | 473 | /* Top left corner */ 474 | if (xints == 2) 475 | { 476 | i = (my-gys-1)*gxm; /* Array index of the top left corner point */ 477 | f[C(i)] = ch_residual_2d 478 | (x, alpha, beta, mparam, hxm2, hym2, 479 | C(i-2*gxm), C(i-gxm+1), C(i-gxm), C(i-gxm+1), 480 | C(i+2), C(i+1), C(i), C(i+1), C(i+2), 481 | C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-2*gxm)); 482 | f[C(i+1)] = ch_residual_2d 483 | (x, alpha, beta, mparam, hxm2, hym2, 484 | C(i-2*gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2), 485 | C(i+1), C(i), C(i+1), C(i+2), C(i+3), 486 | C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-2*gxm+1)); 487 | f[C(i-gxm)] = ch_residual_2d 488 | (x, alpha, beta, mparam, hxm2, hym2, 489 | C(i-gxm), C(i+1), C(i), C(i+1), 490 | C(i-gxm+2), C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2), 491 | C(i-2*gxm+1), C(i-2*gxm), C(i-2*gxm+1), C(i-3*gxm)); 492 | f[C(i-gxm+1)] = ch_residual_2d 493 | (x, alpha, beta, mparam, hxm2, hym2, 494 | C(i-gxm+1), C(i), C(i+1), C(i+2), 495 | C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-gxm+3), 496 | C(i-2*gxm), C(i-2*gxm+1), C(i-2*gxm+2), C(i-3*gxm+1)); 497 | } 498 | 499 | /* Top right corner */ 500 | if (xinte == mx-2) 501 | { /* Array index of top right corner point */ 502 | i = (my-gys-1)*gxm + mx-gxs-1; 503 | f[C(i)] = ch_residual_2d 504 | (x, alpha, beta, mparam, hxm2, hym2, 505 | C(i-2*gxm), C(i-gxm-1), C(i-gxm), C(i-gxm-1), 506 | C(i-2), C(i-1), C(i), C(i-1), C(i-2), 507 | C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-2*gxm)); 508 | f[C(i-1)] = ch_residual_2d 509 | (x, alpha, beta, mparam, hxm2, hym2, 510 | C(i-2*gxm-1), C(i-gxm-2), C(i-gxm-1), C(i-gxm), 511 | C(i-3), C(i-2), C(i-1), C(i), C(i-1), 512 | C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-2*gxm-1)); 513 | f[C(i-gxm)] = ch_residual_2d 514 | (x, alpha, beta, mparam, hxm2, hym2, 515 | C(i-gxm), C(i-1), C(i), C(i-1), 516 | C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-gxm-2), 517 | C(i-2*gxm-1), C(i-2*gxm), C(i-2*gxm-1), C(i-3*gxm)); 518 | f[C(i-gxm-1)] = ch_residual_2d 519 | (x, alpha, beta, mparam, hxm2, hym2, 520 | C(i-gxm-1), C(i-2), C(i-1), C(i), 521 | C(i-gxm-3), C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-gxm-1), 522 | C(i-2*gxm-2), C(i-2*gxm-1), C(i-2*gxm), C(i-3*gxm-1)); 523 | } 524 | } 525 | 526 | /* Test whether we are on the left edge of the global array */ 527 | if (xints == 2) 528 | { 529 | printf ("This must only be reached if we are NOT x-periodic\n"); 530 | /* Left 2 edge lines */ 531 | for (j=yints-gys; j<yinte-gys; j++) 532 | { 533 | i = j*gxm; 534 | f[C(i)] = ch_residual_2d 535 | (x, alpha, beta, mparam, hxm2, hym2, 536 | C(i+2*gxm), C(i+gxm+1), C(i+gxm), C(i+gxm+1), 537 | C(i+2), C(i+1), C(i), C(i+1), C(i+2), 538 | C(i-gxm+1), C(i-gxm), C(i-gxm+1), C(i-2*gxm)); 539 | f[C(i+1)] = ch_residual_2d 540 | (x, alpha, beta, mparam, hxm2, hym2, 541 | C(i+2*gxm+1), C(i+gxm), C(i+gxm+1), C(i+gxm+2), 542 | C(i+1), C(i), C(i+1), C(i+2), C(i+3), 543 | C(i-gxm), C(i-gxm+1), C(i-gxm+2), C(i-2*gxm+1)); 544 | } 545 | } 546 | 547 | /* Test whether we are on the right edge of the global array */ 548 | if (xinte == mx-2) 549 | { 550 | printf ("This must only be reached if we are NOT x-periodic\n"); 551 | /* Right 2 edge lines */ 552 | for (j=yints-gys; j<yinte-gys; j++) 553 | { 554 | i = j*gxm + mx-gxs-1; 555 | f[C(i)] = ch_residual_2d 556 | (x, alpha, beta, mparam, hxm2, hym2, 557 | C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm-1), 558 | C(i-2), C(i-1), C(i), C(i-1), C(i-2), 559 | C(i-gxm-1), C(i-gxm), C(i-gxm-1), C(i-2*gxm)); 560 | f[C(i-1)] = ch_residual_2d 561 | (x, alpha, beta, mparam, hxm2, hym2, 562 | C(i+2*gxm-1), C(i+gxm-2), C(i+gxm-1), C(i+gxm), 563 | C(i-3), C(i-2), C(i-1), C(i), C(i-1), 564 | C(i-gxm-2), C(i-gxm-1), C(i-gxm), C(i-2*gxm-1)); 565 | } 566 | } 567 | 568 | /* Restore vectors */ 569 | ierr = VecRestoreArray (localX, &x); CHKERRQ (ierr); 570 | ierr = VecRestoreArray (localF, &f); CHKERRQ (ierr); 571 | 572 | /* Insert values into global vector */ 573 | ierr = DALocalToGlobal(user->da,localF,INSERT_VALUES,F); CHKERRQ (ierr); 574 | 575 | /* Flop count (multiply-adds are counted as 2 operations) */ 576 | ierr = PetscLogFlops(RESIDUAL_FLOPS_2D*ym*xm); CHKERRQ (ierr); 577 | 578 | return 0; 579 | } 580 | 581 | 582 | /*++++++++++++++++++++++++++++++++++++++ 583 | This evaluates the nonlinear finite difference approximation to the residuals 584 | +latex+$R_i$. 585 | +html+ <i>R<sub>i</sub></i>. 586 | Note that it loops on the interior points and the boundary separately, to 587 | avoid conditional statements within the double loop over the local grid 588 | indices. 589 | +latex+\par 590 | +html+ <p> 591 | Thus far, only periodic boundary conditions work. 592 | 593 | int ch_residual_vector_3d It returns zero or an error value. 594 | 595 | Vec X The pre-allocated local vector of unknowns. 596 | 597 | Vec F The pre-allocated local vector of residuals, filled by this function. 598 | 599 | void *ptr Data passed in the application context. 600 | ++++++++++++++++++++++++++++++++++++++*/ 601 | 602 | int ch_residual_vector_3d (Vec X, Vec F, void *ptr) 603 | { 604 | AppCtx *user = (AppCtx*)ptr; 605 | int ierr,i,j,k, mc,chvar, mx,my,mz,xs,ys,zs,xm,ym,zm; 606 | int gxs,gys,gzs,gxm,gym,gzm, xints,xinte,yints,yinte,zints,zinte; 607 | PetscScalar hx,hy,hz,dhx,dhy,dhz,hxm2,hym2,hzm2; 608 | PetscScalar alpha, beta, mparam; 609 | PetscScalar *x,*f; 610 | Vec localX = user->localX,localF = user->localF; 611 | 612 | /* Scatter ghost points to local vector, using the 2-step process 613 | DAGlobalToLocalBegin(), DAGlobalToLocalEnd(). 614 | By placing code between these two statements, computations can be 615 | done while messages are in transition. */ 616 | ierr = DAGlobalToLocalBegin (user->da, X, INSERT_VALUES, localX); 617 | CHKERRQ (ierr); 618 | 619 | /* Define mesh intervals ratios for uniform grid. 620 | [Note: FD formulae below may someday be normalized by multiplying through 621 | by local volume element to obtain coefficients O(1) in two dimensions.] */ 622 | 623 | mc = user->mc; 624 | chvar = user->chvar; 625 | mx = user->mx; my = user->my; mz = user->mz; 626 | mparam = user->mparam; 627 | alpha = user->kappa * user->gamma * user->epsilon; 628 | beta = user->kappa * user->gamma / user->epsilon; 629 | dhx = (PetscScalar)(mx-1); dhy = (PetscScalar)(my-1); 630 | dhz = (PetscScalar)(mz-1); 631 | /* This next line hard-codes the 1x1x1 cube */ 632 | hx = 1./dhx; hy = 1./dhy; hz = 1./dhz; 633 | hxm2 = 1./hx/hx; hym2 = 1./hy/hy; hzm2 = 1./hz/hz; 634 | 635 | ierr = DAGlobalToLocalEnd (user->da, X, INSERT_VALUES, localX); 636 | CHKERRQ (ierr); 637 | 638 | /* Get pointers to vector data */ 639 | ierr = VecGetArray (localX, &x); CHKERRQ (ierr); 640 | ierr = VecGetArray (localF, &f); CHKERRQ (ierr); 641 | 642 | /* Get local grid boundaries */ 643 | ierr = DAGetCorners (user->da, &xs, &ys, &zs, &xm, &ym, &zm); CHKERRQ (ierr); 644 | ierr = DAGetGhostCorners (user->da, &gxs, &gys, &gzs, &gxm, &gym, &gzm); 645 | CHKERRQ (ierr); 646 | 647 | /* Determine the interior of the locally owned part of the grid. */ 648 | if (user->period == DA_XYZPERIODIC) { 649 | xints = xs; xinte = xs+xm; 650 | yints = ys; yinte = ys+ym; 651 | zints = zs; zinte = zs+zm; } 652 | else { 653 | SETERRQ(PETSC_ERR_ARG_WRONGSTATE, 654 | "Only periodic boundary conditions in 3-D Cahn-Hilliard"); } 655 | 656 | /* Loop over the interior points (no boundaries yet) */ 657 | for (k=zints-gzs; k<zinte-gzs; k++) 658 | for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++) 659 | for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) 660 | f[C(i)] = ch_residual_3d 661 | (x, alpha, beta, mparam, hxm2, hym2, hzm2, 662 | C(i+2*gxm*gym), C(i+gxm*gym+gxm), C(i+gxm*gym-1), C(i+gxm*gym), 663 | C(i+gxm*gym+1), C(i+gxm*gym-gxm), 664 | C(i+2*gxm), C(i+gxm-1), C(i+gxm), C(i+gxm+1), 665 | C(i-2), C(i-1), C(i), C(i+1), C(i+2), 666 | C(i-gxm-1), C(i-gxm), C(i-gxm+1), C(i-2*gxm), 667 | C(i-gxm*gym+gxm), C(i-gxm*gym-1), C(i-gxm*gym), C(i-gxm*gym+1), 668 | C(i-gxm*gym-gxm), C(i-2*gxm*gym)); 669 | 670 | /* Restore vectors */ 671 | ierr = VecRestoreArray (localX, &x); CHKERRQ (ierr); 672 | ierr = VecRestoreArray (localF, &f); CHKERRQ (ierr); 673 | 674 | /* Insert values into global vector */ 675 | ierr = DALocalToGlobal(user->da,localF,INSERT_VALUES,F); CHKERRQ (ierr); 676 | 677 | /* Flop count (multiply-adds are counted as 2 operations) */ 678 | ierr = PetscLogFlops(RESIDUAL_FLOPS_3D*ym*xm*zm); CHKERRQ (ierr); 679 | 680 | /* ierr = VecView (F, VIEWER_STDOUT_SELF); CHKERRQ (ierr); */ 681 | 682 | return 0; 683 | } 684 | 685 | 686 | /*++++++++++++++++++++++++++++++++++++++ 687 | This computes the Jacobian matrix at each iteration, starting with the alpha 688 | term which is pre-computed at the beginning by 689 | +latex+{\tt ch\_jacobian\_alpha\_2d()}. 690 | +html+ <tt>ch_jacobian_alpha_2d()</tt>. 691 | 692 | int ch_jacobian_2d It returns 0 or an error code. 693 | 694 | Vec X The vector of unknowns. 695 | 696 | Mat *A The Jacobian matrix returned to PETSc. 697 | 698 | Mat *B The matrix preconditioner, in this case the same matrix. 699 | 700 | MatStructure *flag Flag saying the nonzeroes are the same each time. 701 | 702 | void *ptr Application context structure. 703 | ++++++++++++++++++++++++++++++++++++++*/ 704 | 705 | int ch_jacobian_2d (Vec X, Mat *A, Mat *B, MatStructure *flag, void *ptr) 706 | { 707 | AppCtx *user = (AppCtx*)ptr; 708 | int ierr,node,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym; 709 | int xints,xinte,yints,yinte; 710 | int columns [5]; 711 | PetscScalar hx,hy,dhx,dhy,hxm2,hym2; 712 | PetscScalar alpha,beta,mparam; 713 | PetscScalar *x, bvalue[5]; 714 | Vec localX = user->localX; 715 | 716 | /* Set the MatStructure flag right, Mats to return */ 717 | *flag = SAME_NONZERO_PATTERN; 718 | A = &(user->J); 719 | B = &(user->J); 720 | 721 | /* Copy the alpha term Jacobian into the main Jacobian space */ 722 | ierr = MatCopy (user->alpha, user->J, SAME_NONZERO_PATTERN); 723 | 724 | /* Scatter ghost points to local vector, using 2-step async I/O with (a 725 | couple of trivial) computations in between. */ 726 | 727 | ierr = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr); 728 | mc = user->mc; chvar = user->chvar; 729 | mx = user->mx; my = user->my; mparam = user->mparam; 730 | alpha = user->kappa * user->gamma * user->epsilon; 731 | beta = user->kappa * user->gamma / user->epsilon; 732 | 733 | /* Define mesh intervals ratios for uniform grid. 734 | [Note: FD formulae below may someday be normalized by multiplying through 735 | by local volume element to obtain coefficients O(1) in two dimensions.] */ 736 | dhx = (PetscScalar)(mx-1); dhy = (PetscScalar)(my-1); 737 | hx = 1./dhx; hy = 1./dhy; /* This line hard-codes the 1x1 square */ 738 | hxm2 = 1./hx/hx; hym2 = 1./hy/hy; 739 | 740 | ierr = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr); 741 | 742 | /* Get pointer to vector data */ 743 | ierr = VecGetArray(localX,&x); CHKERRQ (ierr); 744 | 745 | /* 746 | Get local grid boundaries 747 | */ 748 | ierr = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ (ierr); 749 | ierr = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ (ierr); 750 | 751 | /* Determine the interior of the locally owned part of the grid. */ 752 | if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) { 753 | xints = xs; xinte = xs+xm; } 754 | else { 755 | xints = (xs==0) ? 1:xs; xinte = (xs+xm==mx) ? xs+xm-1 : xs+xm; } 756 | if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) { 757 | yints = ys; yinte = ys+ym; } 758 | else { 759 | yints = (ys==0) ? 1:ys; yinte = (ys+ym==my) ? ys+ym-1 : ys+ym; } 760 | 761 | /* Loop over the interior points... */ 762 | for (j=yints-gys; j<yinte-gys; j++) 763 | for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) { 764 | 765 | /* Form the column indices */ 766 | columns[0]=C(i+gxm); 767 | columns[1]=C(i-1); columns[2]=C(i); columns[3]=C(i+1); 768 | columns[4]=C(i-gxm); 769 | 770 | bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam); 771 | bvalue[1] = beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam); 772 | bvalue[2] = -2.*beta * (hxm2+hym2) * 773 | ch_psidoubleprime (x[columns[2]], mparam); 774 | bvalue[3] = beta * hxm2 * ch_psidoubleprime (x[columns[3]], mparam); 775 | bvalue[4] = beta * hym2 * ch_psidoubleprime (x[columns[4]], mparam); 776 | 777 | node = C(i); 778 | MatSetValuesLocal (user->J, 1,&node, 5,columns, bvalue, ADD_VALUES); 779 | } 780 | 781 | /* Now the boundary conditions... */ 782 | 783 | /* Test whether we're on the bottom row and non-y-periodic */ 784 | if (yints == 1) { 785 | printf ("This must only be reached if we are NOT y-periodic\n"); 786 | for (i=xints-gxs; i<xinte-gxs; i++) { 787 | 788 | /* Form the column indices */ 789 | columns[0]=C(i+gxm); 790 | columns[1]=C(i-1); columns[2]=C(i); columns[3]=C(i+1); 791 | 792 | bvalue[0] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam); 793 | bvalue[1] = beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam); 794 | bvalue[2] = -2.*beta * (hxm2+hym2) * 795 | ch_psidoubleprime (x[columns[2]], mparam); 796 | bvalue[3] = beta * hxm2 * ch_psidoubleprime (x[columns[3]], mparam); 797 | 798 | node = C(i); 799 | MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES); 800 | } 801 | 802 | /* Bottom left corner */ 803 | if (xints == 1) { 804 | columns[0]=C(i=0); columns[1]=C(i+1); 805 | columns[2]=C(i+gxm); 806 | 807 | bvalue[0] = -2.*beta * (hxm2+hym2) * 808 | ch_psidoubleprime (x[columns[0]], mparam); 809 | bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam); 810 | bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam); 811 | 812 | node = C(i); 813 | MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES); 814 | } 815 | 816 | /* Bottom right corner */ 817 | if (xinte == mx-1) { 818 | columns[0]=C(i=mx-gxs-1); columns[1]=C(i-1); 819 | columns[2]=C(i+gxm); 820 | 821 | bvalue[0] = -2.*beta * (hxm2+hym2) * 822 | ch_psidoubleprime (x[columns[0]], mparam); 823 | bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam); 824 | bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam); 825 | 826 | node = C(i); 827 | MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES); 828 | } 829 | } 830 | 831 | /* Test whether we're on the top row and non-y-periodic */ 832 | if (yinte == my-1) { 833 | printf ("This must only be reached if we are NOT y-periodic\n"); 834 | for (i=(my-gys-1)*gxm + xints-gxs; i<(my-gys-1)*gxm + xinte-gxs; i++) { 835 | 836 | /* Form the column indices */ 837 | columns[0]=C(i-1); columns[1]=C(i); columns[2]=C(i+1); 838 | columns[3]=C(i-gxm); 839 | 840 | bvalue[0] = beta * hxm2 * ch_psidoubleprime (x[columns[0]], mparam); 841 | bvalue[1] = -2.*beta * (hxm2+hym2) * 842 | ch_psidoubleprime (x[columns[1]], mparam); 843 | bvalue[2] = beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam); 844 | bvalue[3] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam); 845 | 846 | node = C(i); 847 | MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES); 848 | } 849 | 850 | /* Top left corner */ 851 | if (xints == 1) { 852 | columns[0]=C(i=(my-gys-1)*gxm); columns[1] = C(i+1); 853 | columns[2]=C(i-gxm); 854 | 855 | bvalue[0] = -2.*beta * (hxm2+hym2) * 856 | ch_psidoubleprime (x[columns[0]], mparam); 857 | bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam); 858 | bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam); 859 | 860 | node = C(i); 861 | MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES); 862 | } 863 | 864 | /* Top right corner */ 865 | if (xinte == mx-1) { 866 | columns[0]=C(i=(my-gys-1)*gxm + mx-gxs-1); columns[1]=C(i-1); 867 | columns[2]=C(i-gxm); 868 | 869 | bvalue[0] = -2.*beta * (hxm2+hym2) * 870 | ch_psidoubleprime (x[columns[0]], mparam); 871 | bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam); 872 | bvalue[2] = 2.*beta * hym2 * ch_psidoubleprime (x[columns[2]], mparam); 873 | 874 | node = C(i); 875 | MatSetValuesLocal (user->J, 1,&node, 3,columns, bvalue, ADD_VALUES); 876 | } 877 | } 878 | 879 | /* Left column */ 880 | if (xints == 1) { 881 | printf ("This must only be reached if we are NOT x-periodic\n"); 882 | for (i=(yints-gys)*gxm; i<(yinte-gys)*gxm; i+=gxm) { 883 | columns[0]=C(i+gxm); 884 | columns[1]=C(i); columns[2] = C(i+1); 885 | columns[3]=C(i-gxm); 886 | 887 | bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam); 888 | bvalue[1] = -2.*beta * (hxm2+hym2) * 889 | ch_psidoubleprime (x[columns[1]], mparam); 890 | bvalue[2] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam); 891 | bvalue[3] = beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam); 892 | 893 | node = C(i); 894 | MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES); 895 | } 896 | } 897 | 898 | /* Right column */ 899 | if (xinte == mx-1) { 900 | printf ("This must only be reached if we are NOT x-periodic\n"); 901 | for (i=(yints-gys)*gxm + mx-gxs-1; i<(yinte-gys)*gxm + mx-gxs-1; i+=gxm) { 902 | columns[0]=C(i+gxm); 903 | columns[1]=C(i-1); columns[2] = C(i); 904 | columns[3]=C(i-gxm); 905 | 906 | bvalue[0] = beta * hym2 * ch_psidoubleprime (x[columns[0]], mparam); 907 | bvalue[1] = 2.*beta * hxm2 * ch_psidoubleprime (x[columns[1]], mparam); 908 | bvalue[2] = -2.*beta * (hxm2+hym2) * 909 | ch_psidoubleprime (x[columns[2]], mparam); 910 | bvalue[3] = beta * hym2 * ch_psidoubleprime (x[columns[3]], mparam); 911 | 912 | node = C(i); 913 | MatSetValuesLocal (user->J, 1,&node, 4,columns, bvalue, ADD_VALUES); 914 | } 915 | } 916 | 917 | /* Assemble matrix, using the 2-step process: 918 | MatAssemblyBegin(), MatAssemblyEnd(). */ 919 | ierr = MatAssemblyBegin(user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 920 | ierr = MatAssemblyEnd(user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 921 | 922 | /* Restore unknown vector */ 923 | ierr = VecRestoreArray(localX,&x); CHKERRQ (ierr); 924 | 925 | return ierr; 926 | } 927 | 928 | 929 | /*++++++++++++++++++++++++++++++++++++++ 930 | This computes the Jacobian matrix at each iteration, starting with the alpha 931 | term which is pre-computed at the beginning by 932 | +latex+{\tt ch\_jacobian\_alpha\_3d()}. 933 | +html+ <tt>ch_jacobian_alpha_3d()</tt>. 934 | 935 | int ch_jacobian_3d It returns 0 or an error code. 936 | 937 | Vec X The vector of unknowns. 938 | 939 | Mat *A The Jacobian matrix returned to PETSc. 940 | 941 | Mat *B The matrix preconditioner, in this case the same matrix. 942 | 943 | MatStructure *flag Flag saying the nonzeroes are the same each time. 944 | 945 | void *ptr Application context structure. 946 | ++++++++++++++++++++++++++++++++++++++*/ 947 | 948 | int ch_jacobian_3d (Vec X, Mat *A, Mat *B, MatStructure *flag, void *ptr) 949 | { 950 | AppCtx *user = (AppCtx*)ptr; 951 | int ierr,node,i,j,k, mc,chvar, mx,my,mz, xs,ys,zs, xm,ym,zm; 952 | int gxs,gys,gzs, gxm,gym,gzm, xints,xinte, yints,yinte, zints,zinte; 953 | int columns [7]; 954 | PetscScalar hx,hy,hz, dhx,dhy,dhz, hxm2,hym2,hzm2; 955 | PetscScalar alpha,beta,mparam; 956 | PetscScalar *x, bvalue[7]; 957 | Vec localX = user->localX; 958 | 959 | /* Set the MatStructure flag right, Mats to return */ 960 | *flag = SAME_NONZERO_PATTERN; 961 | A = &(user->J); 962 | B = &(user->J); 963 | 964 | /* Copy the alpha term Jacobian into the main Jacobian space */ 965 | ierr = MatCopy (user->alpha, user->J, SAME_NONZERO_PATTERN); 966 | 967 | /* Scatter ghost points to local vector, using 2-step async I/O with (a 968 | couple of trivial) computations in between. */ 969 | 970 | ierr = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr); 971 | mc = user->mc; chvar = user->chvar; 972 | mx = user->mx; my = user->my; mz = user->mz; mparam = user->mparam; 973 | alpha = user->kappa * user->gamma * user->epsilon; 974 | beta = user->kappa * user->gamma / user->epsilon; 975 | 976 | /* Define mesh intervals ratios for uniform grid. 977 | [Note: FD formulae below may someday be normalized by multiplying through 978 | by local volume element to obtain coefficients O(1) in two dimensions.] */ 979 | dhx = (PetscScalar)(mx-1); dhy = (PetscScalar)(my-1); 980 | dhz = (PetscScalar)(mz-1); 981 | /* This line hard-codes the 1x1x1 cube */ 982 | hx = 1./dhx; hy = 1./dhy; hz = 1./dhz; 983 | hxm2 = 1./hx/hx; hym2 = 1./hy/hy; hzm2 = 1./hz/hz; 984 | 985 | ierr = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ (ierr); 986 | 987 | /* Get pointer to vector data */ 988 | ierr = VecGetArray(localX,&x); CHKERRQ (ierr); 989 | 990 | /* Get local grid boundaries */ 991 | ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr); 992 | ierr = DAGetGhostCorners(user->da,&gxs,&gys,&gzs,&gxm,&gym,&gzm); 993 | CHKERRQ (ierr); 994 | 995 | /* Determine the interior of the locally owned part of the grid. */ 996 | if (user->period == DA_XYZPERIODIC) { 997 | xints = xs; xinte = xs+xm; 998 | yints = ys; yinte = ys+ym; 999 | zints = zs; zinte = zs+zm; } 1000 | else { 1001 | SETERRQ(PETSC_ERR_ARG_WRONGSTATE, 1002 | "Only periodic boundary conditions in 3-D Cahn-Hilliard"); } 1003 | 1004 | /* Loop over the interior points... */ 1005 | for (k=zints-gzs; k<zinte-gzs; k++) 1006 | for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++) 1007 | for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) { 1008 | 1009 | /* Form the column indices */ 1010 | columns[0]=C(i+gxm*gym); columns[1]=C(i+gxm); 1011 | columns[2]=C(i-1); columns[3]=C(i); columns[4]=C(i+1); 1012 | columns[5]=C(i-gxm); columns[6]=C(i-gxm*gym); 1013 | 1014 | bvalue[0] = beta * hzm2 * ch_psidoubleprime (x[columns[0]], mparam); 1015 | bvalue[1] = beta * hym2 * ch_psidoubleprime (x[columns[1]], mparam); 1016 | bvalue[2] = beta * hxm2 * ch_psidoubleprime (x[columns[2]], mparam); 1017 | bvalue[3] = -2.*beta * (hxm2+hym2+hzm2) * 1018 | ch_psidoubleprime (x[columns[3]], mparam); 1019 | bvalue[4] = beta * hxm2 * ch_psidoubleprime (x[columns[4]], mparam); 1020 | bvalue[5] = beta * hym2 * ch_psidoubleprime (x[columns[5]], mparam); 1021 | bvalue[6] = beta * hzm2 * ch_psidoubleprime (x[columns[6]], mparam); 1022 | 1023 | node = C(i); 1024 | MatSetValuesLocal (user->J, 1,&node, 5,columns, bvalue, ADD_VALUES); 1025 | } 1026 | 1027 | /* Assemble matrix, using the 2-step process: 1028 | MatAssemblyBegin(), MatAssemblyEnd(). */ 1029 | ierr = MatAssemblyBegin (user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 1030 | ierr = MatAssemblyEnd (user->J,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 1031 | 1032 | /* Restore unknown vector */ 1033 | ierr = VecRestoreArray (localX,&x); CHKERRQ (ierr); 1034 | 1035 | /* ierr = MatView (user->J, VIEWER_STDOUT_SELF); CHKERRQ (ierr); */ 1036 | 1037 | return ierr; 1038 | } 1039 | 1040 | 1041 | /*++++++++++++++++++++++++++++++++++++++ 1042 | This creates the initial alpha and J matrices, where alpha is the alpha term 1043 | component of the Jacobian. Since the alpha term is linear, this part of the 1044 | Jacobian need only be calculated once. 1045 | 1046 | int ch_jacobian_alpha_2d It returns zero or an error message. 1047 | 1048 | AppCtx *user The application context structure pointer. 1049 | ++++++++++++++++++++++++++++++++++++++*/ 1050 | 1051 | int ch_jacobian_alpha_2d (AppCtx *user) 1052 | { 1053 | int ierr,node,i,j, mc,chvar, mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym; 1054 | int xints,xinte,yints,yinte; 1055 | int columns [13]; 1056 | PetscScalar hx,hy,dhx,dhy,hxm2,hym2; 1057 | PetscScalar alpha,beta,mparam; 1058 | PetscScalar avalue [25]; 1059 | 1060 | mc = user->mc; chvar = user->chvar; 1061 | mx = user->mx; my = user->my; mparam = user->mparam; 1062 | alpha = user->kappa * user->gamma * user->epsilon; 1063 | beta = user->kappa * user->gamma / user->epsilon; 1064 | 1065 | /* Define mesh intervals ratios for uniform grid. 1066 | [Note: FD formulae below may someday be normalized by multiplying through 1067 | by local volume element to obtain coefficients O(1) in two dimensions.] */ 1068 | dhx = (PetscScalar)(mx-1); dhy = (PetscScalar)(my-1); 1069 | hx = 1./dhx; hy = 1./dhy; /* This line hard-codes the 1x1 square */ 1070 | hxm2 = 1./hx/hx; hym2 = 1./hy/hy; 1071 | 1072 | /* Get local grid boundaries */ 1073 | ierr = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 1074 | CHKERRQ (ierr); 1075 | ierr = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 1076 | CHKERRQ (ierr); 1077 | 1078 | /* Determine the interior of the locally owned part of the grid. */ 1079 | if (user->period == DA_XPERIODIC || user->period == DA_XYPERIODIC) { 1080 | xints = xs; xinte = xs+xm; } 1081 | else { 1082 | xints = (xs==0) ? 2:xs; xinte = (xs+xm==mx) ? xs+xm-2 : xs+xm; } 1083 | if (user->period == DA_YPERIODIC || user->period == DA_XYPERIODIC) { 1084 | yints = ys; yinte = ys+ym; } 1085 | else { 1086 | yints = (ys==0) ? 2:ys; yinte = (ys+ym==my) ? ys+ym-2 : ys+ym; } 1087 | 1088 | /* Get parameters for matrix creation */ 1089 | i = xm * ym * user->mc; 1090 | j = mx * my * user->mc; 1091 | 1092 | /* Create the distributed alpha matrix */ 1093 | ierr = MatCreateMPIAIJ (user->comm, i,i, j,j, 13,PETSC_NULL, 6,PETSC_NULL, 1094 | &(user->alpha)); 1095 | CHKERRQ (ierr); 1096 | 1097 | /* Get the local-to-global mapping and associate it with the matrix */ 1098 | ierr = DAGetISLocalToGlobalMapping (user->da, &user->isltog); CHKERRQ (ierr); 1099 | ierr = MatSetLocalToGlobalMapping (user->alpha, user->isltog); CHKERRQ (ierr); 1100 | 1101 | /* Pre-compute alpha term values (see ch_residual_2d above) 1102 | This should be the only place they're actually computed; they'll be 1103 | used below. */ 1104 | avalue[0] = avalue[12] = -alpha*hym2*hym2; 1105 | avalue[1] = avalue[3] = avalue[9] = avalue[11] = -2.*alpha*hxm2*hym2; 1106 | avalue[2] = avalue[10] = 4.*alpha*hym2*(hxm2+hym2); 1107 | avalue[4] = avalue[8] = -alpha*hxm2*hxm2; 1108 | avalue[5] = avalue[7] = 4.*alpha*hxm2*(hxm2+hym2); 1109 | avalue[6] = -alpha*(6.*hxm2*hxm2 + 6.*hym2*hym2 + 8.*hxm2*hym2); 1110 | 1111 | /* Loop over interior nodes */ 1112 | for (j=yints-gys; j<yinte-gys; j++) 1113 | for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) { 1114 | 1115 | /* Form the column indices */ 1116 | columns[0]=C(i+2*gxm); 1117 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1118 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1119 | columns[7]=C(i+1); columns[8]=C(i+2); 1120 | columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1); 1121 | columns[12]=C(i-2*gxm); 1122 | 1123 | node = C(i); 1124 | MatSetValuesLocal (user->alpha, 1,&node, 13,columns, avalue, 1125 | INSERT_VALUES); 1126 | } 1127 | 1128 | /* Okay, that was the easy part, now for the boundary conditions! */ 1129 | 1130 | /* Test whether we're on the bottom row and non-y-periodic */ 1131 | if (yints == 2) { 1132 | printf ("This must only be reached if we are NOT y-periodic\n"); 1133 | /* Change value 6 and do the second-from-bottom row */ 1134 | avalue[6] += avalue[12]; 1135 | for (i=gxm + xints-gxs; i<gxm + xinte-gxs; i++) { 1136 | 1137 | /* Form the column indices */ 1138 | columns[0]=C(i+2*gxm); 1139 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1140 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1141 | columns[7]=C(i+1); columns[8]=C(i+2); 1142 | columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1); 1143 | 1144 | node = C(i); 1145 | MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue, 1146 | INSERT_VALUES); 1147 | } 1148 | avalue[6] -= avalue[12]; 1149 | 1150 | /* Make some new avalues and do the bottom row */ 1151 | avalue[13] = 2.*avalue[0]; 1152 | avalue[14] = avalue[16] = 2.*avalue[1]; 1153 | avalue[15] = 2.*avalue[2]; 1154 | avalue[17] = avalue[21] = avalue[4]; 1155 | avalue[18] = avalue[20] = avalue[5]; 1156 | avalue[19] = avalue[6]; 1157 | for (i=xints-gxs; i<xinte-gxs; i++) { 1158 | 1159 | /* Form the column indices */ 1160 | columns[0]=C(i+2*gxm); 1161 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1162 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1163 | columns[7]=C(i+1); columns[8]=C(i+2); 1164 | 1165 | node = C(i); 1166 | MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13, 1167 | INSERT_VALUES); 1168 | } 1169 | 1170 | /* Bottom left corner */ 1171 | if (xints == 2) { 1172 | node=C(0); /* The point (0,0) */ 1173 | columns[0]=C(2*gxm); 1174 | columns[1]=C(gxm); columns[2]=C(gxm+1); 1175 | columns[3]=C(0); columns[4]=C(1); columns[5]=C(2); 1176 | avalue[13] = 2.*avalue[0]; 1177 | avalue[14] = 2.*avalue[2]; 1178 | avalue[15] = 4.*avalue[3]; 1179 | avalue[16] = avalue[6]; 1180 | avalue[17] = 2.*avalue[7]; 1181 | avalue[18] = 2.*avalue[8]; 1182 | MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13, 1183 | INSERT_VALUES); 1184 | node=C(1); /* The point (1,0) */ 1185 | columns[0]=C(2*gxm+1); 1186 | columns[1]=C(gxm); columns[2]=C(gxm+1); columns[3]=C(gxm+2); 1187 | columns[4]=C(0); columns[5]=C(1); columns[6]=C(2); columns[7]=C(3); 1188 | avalue[13] = 2.*avalue[0]; 1189 | avalue[14] = avalue[16] = 2.*avalue[1]; 1190 | avalue[15] = 2.*avalue[2]; 1191 | avalue[17] = avalue[19] = avalue[5]; 1192 | avalue[18] = avalue[6] + avalue[4]; 1193 | avalue[20] = avalue[8]; 1194 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1195 | INSERT_VALUES); 1196 | node=C(gxm); /* The point (0,1) */ 1197 | columns[0]=C(3*gxm); 1198 | columns[1]=C(2*gxm); columns[2]=C(2*gxm+1); 1199 | columns[3]=C(gxm); columns[4]=C(gxm+1); columns[5]=C(gxm+2); 1200 | columns[6]=C(0); columns[7]=C(1); 1201 | avalue[13] = avalue[0]; 1202 | avalue[14] = avalue[19] = avalue[2]; 1203 | avalue[15] = avalue[20] = 2.*avalue[3]; 1204 | avalue[16] = avalue[6] + avalue[12]; 1205 | avalue[17] = 2.*avalue[7]; 1206 | avalue[18] = 2.*avalue[8]; 1207 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1208 | INSERT_VALUES); 1209 | node=C(gxm+1); /* The point (1,1) */ 1210 | columns[0]=C(3*gxm+1); 1211 | columns[1]=C(2*gxm); columns[2]=C(2*gxm+1); columns[3]=C(2*gxm+2); 1212 | columns[4]=C(gxm); columns[5]=C(gxm+1); columns[6]=C(gxm+2); 1213 | columns[7]=C(gxm+3); 1214 | columns[8]=C(0); columns[9]=C(1); columns[10]=C(2); 1215 | avalue[13] = avalue[0]; 1216 | avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1]; 1217 | avalue[15] = avalue[22] = avalue[2]; 1218 | avalue[17] = avalue[19] = avalue[5]; 1219 | avalue[18] = avalue[6] + avalue[4] + avalue[12]; 1220 | avalue[20] = avalue[8]; 1221 | MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13, 1222 | INSERT_VALUES); 1223 | } 1224 | 1225 | /* Bottom right corner */ 1226 | if(xinte == mx-2) { 1227 | node=C(i=mx-gxs-1); /* The point (mx-1, 0) */ 1228 | columns[0]=C(i+2*gxm); 1229 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); 1230 | columns[3]=C(i-2); columns[4]=C(i-1); columns[5]=C(i); 1231 | avalue[13] = 2.*avalue[0]; 1232 | avalue[14] = 4.*avalue[1]; 1233 | avalue[15] = 2.*avalue[2]; 1234 | avalue[16] = 2.*avalue[4]; 1235 | avalue[17] = 2.*avalue[5]; 1236 | avalue[18] = avalue[6]; 1237 | MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13, 1238 | INSERT_VALUES); 1239 | node=C(i=mx-gxs-2); /* The point (mx-2,0) */ 1240 | columns[0]=C(i+2*gxm); 1241 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1242 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1243 | columns[7]=C(i+1); 1244 | avalue[13] = 2.*avalue[0]; 1245 | avalue[14] = avalue[16] = 2.*avalue[1]; 1246 | avalue[15] = 2.*avalue[2]; 1247 | avalue[17] = avalue[4]; 1248 | avalue[18] = avalue[20] = avalue[5]; 1249 | avalue[19] = avalue[6] + avalue[8]; 1250 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1251 | INSERT_VALUES); 1252 | node=C(i=gxm+mx-gxs-1); /* The point (mx-1,1) */ 1253 | columns[0]=C(i+2*gxm); 1254 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); 1255 | columns[3]=C(i-2); columns[4]=C(i-1); columns[5]=C(i); 1256 | columns[6]=C(i-gxm-1); columns[7]=C(i-gxm); 1257 | avalue[13] = avalue[0]; 1258 | avalue[14] = avalue[19] = 2.*avalue[1]; 1259 | avalue[15] = avalue[20] = avalue[2]; 1260 | avalue[16] = 2.*avalue[4]; 1261 | avalue[17] = 2.*avalue[5]; 1262 | avalue[18] = avalue[6] + avalue[12]; 1263 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1264 | INSERT_VALUES); 1265 | node=C(i=gxm+mx-gxs-2); /* The point (mx-2,1) */ 1266 | columns[0]=C(i+2*gxm); 1267 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1268 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1269 | columns[7]=C(i+1); 1270 | columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1); 1271 | avalue[13] = avalue[0]; 1272 | avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1]; 1273 | avalue[15] = avalue[22] = avalue[2]; 1274 | avalue[17] = avalue[4]; 1275 | avalue[18] = avalue[20] = avalue[5]; 1276 | avalue[19] = avalue[6] + avalue[8] + avalue[12]; 1277 | MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13, 1278 | INSERT_VALUES); 1279 | } 1280 | } 1281 | 1282 | /* Test whether we're on the top row and non-y-periodic */ 1283 | if (yinte == my-2) { 1284 | printf ("This must only be reached if we are NOT y-periodic\n"); 1285 | /* Change value 6 and do the second-from-top row */ 1286 | avalue[6] += avalue[0]; 1287 | for (i=(my-gys-2)*gxm + xints-gxs; 1288 | i<(my-gys-2)*gxm + xinte-gxs; i++) { 1289 | 1290 | /* Form the column indices */ 1291 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1292 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1293 | columns[7]=C(i+1); columns[8]=C(i+2); 1294 | columns[9]=C(i-gxm-1); columns[10]=C(i-gxm); columns[11]=C(i-gxm+1); 1295 | columns[12]=C(i-2*gxm); 1296 | 1297 | node = C(i); 1298 | MatSetValuesLocal (user->alpha, 1,&node, 12,columns+1, avalue+1, 1299 | INSERT_VALUES); 1300 | } 1301 | avalue[6] -= avalue[0]; 1302 | 1303 | /* Make some new avalues and do the top row */ 1304 | avalue[13] = avalue[17] = avalue[4]; 1305 | avalue[14] = avalue[16] = avalue[5]; 1306 | avalue[15] = avalue[6]; 1307 | avalue[18] = avalue[20] = 2.*avalue[9]; 1308 | avalue[19] = 2.*avalue[10]; 1309 | avalue[21] = 2.*avalue[12]; 1310 | for (i=(my-gys-1)*gxm + xints-gxs; 1311 | i<(my-gys-1)*gxm + xinte-gxs; i++) { 1312 | 1313 | /* Form the column indices */ 1314 | columns[0]=C(i-2); columns[1]=C(i-1); columns[2]=C(i); 1315 | columns[3]=C(i+1); columns[4]=C(i+2); 1316 | columns[5]=C(i-gxm-1); columns[6]=C(i-gxm); columns[7]=C(i-gxm+1); 1317 | columns[8]=C(i-2*gxm); 1318 | 1319 | node = C(i); 1320 | MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13, 1321 | INSERT_VALUES); 1322 | } 1323 | 1324 | /* Top left corner */ 1325 | if (xints == 2) { 1326 | node=C(i=(my-gys-1)*gxm); /* The point (0,my-1) */ 1327 | columns[0]=C(i); columns[1]=C(i+1); columns[2]=C(i+2); 1328 | columns[3]=C(i-gxm); columns[4]=C(i-gxm+1); 1329 | columns[5]=C(i-2*gxm); 1330 | avalue[13] = avalue[6]; 1331 | avalue[14] = 2.*avalue[7]; 1332 | avalue[15] = 2.*avalue[8]; 1333 | avalue[16] = 2.*avalue[10]; 1334 | avalue[17] = 4.*avalue[11]; 1335 | avalue[18] = 2.*avalue[12]; 1336 | MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13, 1337 | INSERT_VALUES); 1338 | node=C(i=(my-gys-1)*gxm+1); /* The point (1,my-1) */ 1339 | columns[0]=C(i-1); columns[1]=C(i); columns[2]=C(i+1); 1340 | columns[3]=C(i+2); 1341 | columns[4]=C(i-gxm-1); columns[5]=C(i-gxm); columns[6]=C(i-gxm+1); 1342 | columns[7]=C(i-2*gxm); 1343 | avalue[13] = avalue[15] = avalue[5]; 1344 | avalue[14] = avalue[6] + avalue[4]; 1345 | avalue[16] = avalue[8]; 1346 | avalue[17] = avalue[19] = 2.*avalue[9]; 1347 | avalue[18] = 2.*avalue[10]; 1348 | avalue[20] = 2.*avalue[12]; 1349 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1350 | INSERT_VALUES); 1351 | node=C(i=(my-gys-2)*gxm); /* The point (0,my-2) */ 1352 | columns[0]=C(i+gxm); columns[1]=C(i+gxm+1); 1353 | columns[2]=C(i); columns[3]=C(i+1); columns[4]=C(i+2); 1354 | columns[5]=C(i-gxm); columns[6]=C(i-gxm+1); 1355 | columns[7]=C(i-2*gxm); 1356 | avalue[13] = avalue[18] = avalue[2]; 1357 | avalue[14] = avalue[19] = 2.*avalue[3]; 1358 | avalue[15] = avalue[6] + avalue[0]; 1359 | avalue[16] = 2.*avalue[7]; 1360 | avalue[17] = 2.*avalue[8]; 1361 | avalue[20] = avalue[12]; 1362 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1363 | INSERT_VALUES); 1364 | node=C(i=(my-gys-2)*gxm+1); /* The point (1,my-2) */ 1365 | columns[0]=C(i+gxm-1); columns[1]=C(i+gxm); columns[2]=C(i+gxm+1); 1366 | columns[3]=C(i-1); columns[4]=C(i); columns[5]=C(i+1); 1367 | columns[6]=C(i+2); 1368 | columns[7]=C(i-gxm-1); columns[8]=C(i-gxm); columns[9]=C(i-gxm+1); 1369 | columns[10]=C(i-2*gxm); 1370 | avalue[13] = avalue[15] = avalue[20] = avalue[22] = avalue[1]; 1371 | avalue[14] = avalue[21] = avalue[2]; 1372 | avalue[16] = avalue[18] = avalue[5]; 1373 | avalue[17] = avalue[6] + avalue[4] + avalue[0]; 1374 | avalue[19] = avalue[8]; 1375 | avalue[23] = avalue[12]; 1376 | MatSetValuesLocal (user->alpha, 1,&node, 11,columns, avalue+13, 1377 | INSERT_VALUES); 1378 | } 1379 | 1380 | /* Top right corner */ 1381 | if(xinte == mx-2) { 1382 | node=C(i=(my-gys-1)*gxm + mx-gxs-1); /* The point (mx-1, my-1) */ 1383 | columns[0]=C(i-2); columns[1]=C(i-1); columns[2]=C(i); 1384 | columns[3]=C(i-gxm-1); columns[4]=C(i-gxm); 1385 | columns[5]=C(i-2*gxm); 1386 | avalue[13] = 2.*avalue[4]; 1387 | avalue[14] = 2.*avalue[5]; 1388 | avalue[15] = avalue[6]; 1389 | avalue[16] = 4.*avalue[9]; 1390 | avalue[17] = 2.*avalue[10]; 1391 | avalue[18] = 2.*avalue[12]; 1392 | MatSetValuesLocal (user->alpha, 1,&node, 6,columns, avalue+13, 1393 | INSERT_VALUES); 1394 | node=C(i=(my-gys-1)*gxm + mx-gxs-2); /* The point (mx-2, my-1) */ 1395 | columns[0]=C(i-2); columns[1]=C(i-1); columns[2]=C(i); 1396 | columns[3]=C(i+1); 1397 | columns[4]=C(i-gxm-1); columns[5]=C(i-gxm); columns[6]=C(i-gxm+1); 1398 | columns[7]=C(i-2*gxm); 1399 | avalue[13] = avalue[4]; 1400 | avalue[14] = avalue[16] = avalue[5]; 1401 | avalue[15] = avalue[6] + avalue[8]; 1402 | avalue[17] = avalue[19] = 2.*avalue[9]; 1403 | avalue[18] = 2.*avalue[10]; 1404 | avalue[20] = 2.*avalue[12]; 1405 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1406 | INSERT_VALUES); 1407 | node=C(i=(my-gys-2)*gxm + mx-gxs-1); /* The point (mx-1, my-2) */ 1408 | columns[0]=C(i+gxm-1); columns[1]=C(i+gxm); 1409 | columns[2]=C(i-2); columns[3]=C(i-1); columns[4]=C(i); 1410 | columns[5]=C(i-gxm-1); columns[6]=C(i-gxm); 1411 | columns[7]=C(i-2*gxm); 1412 | avalue[13] = avalue[18] = 2.*avalue[1]; 1413 | avalue[14] = avalue[19] = avalue[2]; 1414 | avalue[15] = 2.*avalue[4]; 1415 | avalue[16] = 2.*avalue[5]; 1416 | avalue[17] = avalue[6] + avalue[0]; 1417 | avalue[20] = avalue[12]; 1418 | MatSetValuesLocal (user->alpha, 1,&node, 8,columns, avalue+13, 1419 | INSERT_VALUES); 1420 | node=C(i=(my-gys-2)*gxm + mx-gxs-2); /* The point (mx-2, my-2) */ 1421 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1422 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1423 | columns[7]=C(i+1); 1424 | columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1); 1425 | columns[11]=C(i-2*gxm); 1426 | avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1]; 1427 | avalue[15] = avalue[22] = avalue[2]; 1428 | avalue[17] = avalue[4]; 1429 | avalue[18] = avalue[20] = avalue[5]; 1430 | avalue[19] = avalue[6] + avalue[8] + avalue[0]; 1431 | avalue[24] = avalue[12]; 1432 | MatSetValuesLocal (user->alpha, 1,&node, 11,columns+1, avalue+14, 1433 | INSERT_VALUES); 1434 | } 1435 | } 1436 | 1437 | /* Test whether we're on the left column and non-x-periodic */ 1438 | if (xints == 2) { 1439 | printf ("This must only be reached if we are NOT x-periodic\n"); 1440 | /* Make some avalues and do the second-from-left column */ 1441 | avalue[13] = avalue[24] = avalue[0]; 1442 | avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1]; 1443 | avalue[15] = avalue[22] = avalue[2]; 1444 | avalue[17] = avalue[19] = avalue[5]; 1445 | avalue[18] = avalue[6] + avalue[4]; 1446 | avalue[20] = avalue[8]; 1447 | for (i=(yints-gys)*gxm + 1; i<(yinte-gys)*gxm + 1; i+=gxm) { 1448 | columns[0]=C(i+2*gxm); 1449 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1450 | columns[4]=C(i-1); columns[5]=C(i); columns[6]=C(i+1); 1451 | columns[7]=C(i+2); 1452 | columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1); 1453 | columns[11]=C(i-2*gxm); 1454 | 1455 | node = C(i); 1456 | MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue+13, 1457 | INSERT_VALUES); 1458 | } 1459 | 1460 | /* Make some more avalues and do the leftmost column */ 1461 | avalue[13] = avalue[21] = avalue[0]; 1462 | avalue[14] = avalue[19] = avalue[2]; 1463 | avalue[15] = avalue[20] = 2.*avalue[3]; 1464 | avalue[16] = avalue[6]; 1465 | avalue[17] = 2.*avalue[7]; 1466 | avalue[18] = 2.*avalue[8]; 1467 | for (i=(yints-gys)*gxm; i<(yinte-gys)*gxm; i+=gxm) { 1468 | columns[0]=C(i+2*gxm); 1469 | columns[1]=C(i+gxm); columns[2]=C(i+gxm+1); 1470 | columns[3]=C(i); columns[4]=C(i+1); columns[5]=C(i+2); 1471 | columns[6]=C(i-gxm); columns[7]=C(i-gxm+1); 1472 | columns[8]=C(i-2*gxm); 1473 | 1474 | node = C(i); 1475 | MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13, 1476 | INSERT_VALUES); 1477 | } 1478 | } 1479 | 1480 | /* Test whether we're on the right column and non-x-periodic */ 1481 | if (xinte == mx-2) { 1482 | printf ("This must only be reached if we are NOT x-periodic\n"); 1483 | /* Make some avalues and do the second-from-right column */ 1484 | avalue[13] = avalue[24] = avalue[0]; 1485 | avalue[14] = avalue[16] = avalue[21] = avalue[23] = avalue[1]; 1486 | avalue[15] = avalue[22] = avalue[2]; 1487 | avalue[17] = avalue[4]; 1488 | avalue[18] = avalue[20] = avalue[5]; 1489 | avalue[19] = avalue[6] + avalue[8]; 1490 | for (i=(yints-gys)*gxm + mx-gxs-2; i<(yinte-gys)*gxm + mx-gxs-2; i+=gxm) { 1491 | columns[0]=C(i+2*gxm); 1492 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); columns[3]=C(i+gxm+1); 1493 | columns[4]=C(i-2); columns[5]=C(i-1); columns[6]=C(i); 1494 | columns[7]=C(i+1); 1495 | columns[8]=C(i-gxm-1); columns[9]=C(i-gxm); columns[10]=C(i-gxm+1); 1496 | columns[11]=C(i-2*gxm); 1497 | 1498 | node = C(i); 1499 | MatSetValuesLocal (user->alpha, 1,&node, 12,columns, avalue+13, 1500 | INSERT_VALUES); 1501 | } 1502 | 1503 | /* Make some more avalues and do the rightmost column */ 1504 | avalue[13] = avalue[21] = avalue[0]; 1505 | avalue[14] = avalue[19] = 2.*avalue[1]; 1506 | avalue[15] = avalue[20] = avalue[2]; 1507 | avalue[16] = 2.*avalue[4]; 1508 | avalue[17] = 2.*avalue[5]; 1509 | avalue[18] = avalue[6]; 1510 | for (i=(yints-gys)*gxm + mx-gxs-1; i<(yinte-gys)*gxm + mx-gxs-1; i+=gxm) { 1511 | columns[0]=C(i+2*gxm); 1512 | columns[1]=C(i+gxm-1); columns[2]=C(i+gxm); 1513 | columns[3]=C(i-2); columns[4]=C(i-1); columns[5]=C(i); 1514 | columns[6]=C(i-gxm-1); columns[7]=C(i-gxm); 1515 | columns[8]=C(i-2*gxm); 1516 | 1517 | node = C(i); 1518 | MatSetValuesLocal (user->alpha, 1,&node, 9,columns, avalue+13, 1519 | INSERT_VALUES); 1520 | } 1521 | } 1522 | 1523 | /* Assemble matrix, using the 2-step process: 1524 | MatAssemblyBegin(), MatAssemblyEnd(). */ 1525 | ierr = MatAssemblyBegin(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 1526 | ierr = MatAssemblyEnd(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 1527 | 1528 | /* Make J a copy of alpha with the same local mapping (setting new mapping is 1529 | unnecessary with PETSc 2.1.5). */ 1530 | ierr = MatDuplicate (user->alpha, MAT_COPY_VALUES, &(user->J)); CHKERRQ (ierr); 1531 | 1532 | /* Tell the matrix we will never add a new nonzero location to the 1533 | matrix. If we do, it will generate an error. */ 1534 | ierr = MatSetOption(user->J,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ (ierr); 1535 | 1536 | return ierr; 1537 | } 1538 | 1539 | 1540 | /*++++++++++++++++++++++++++++++++++++++ 1541 | This creates the initial alpha and J matrices, where alpha is the alpha term 1542 | component of the Jacobian. Since the alpha term is linear, this part of the 1543 | Jacobian need only be calculated once. 1544 | 1545 | int ch_jacobian_alpha_3d It returns zero or an error message. 1546 | 1547 | AppCtx *user The application context structure pointer. 1548 | ++++++++++++++++++++++++++++++++++++++*/ 1549 | 1550 | int ch_jacobian_alpha_3d (AppCtx *user) 1551 | { 1552 | int ierr,node,i,j,k, mc,chvar, mx,my,mz, xs,ys,zs, xm,ym,zm; 1553 | int gxs,gys,gzs, gxm,gym,gzm, xints,xinte, yints,yinte, zints,zinte; 1554 | int columns [25]; 1555 | PetscScalar hx,hy,hz, dhx,dhy,dhz, hxm2,hym2,hzm2; 1556 | PetscScalar alpha,beta,mparam; 1557 | PetscScalar avalue [25]; 1558 | 1559 | mc = user->mc; chvar = user->chvar; 1560 | mx = user->mx; my = user->my; mz = user->mz; mparam = user->mparam; 1561 | alpha = user->kappa * user->gamma * user->epsilon; 1562 | beta = user->kappa * user->gamma / user->epsilon; 1563 | 1564 | /* Define mesh intervals ratios for uniform grid. 1565 | [Note: FD formulae below may someday be normalized by multiplying through 1566 | by local volume element to obtain coefficients O(1) in two dimensions.] */ 1567 | dhx = (PetscScalar)(mx-1); dhy = (PetscScalar)(my-1); 1568 | dhz = (PetscScalar)(mz-1); 1569 | /* This line hard-codes the 1x1x1 cube */ 1570 | hx = 1./dhx; hy = 1./dhy; hz = 1./dhz; 1571 | hxm2 = 1./hx/hx; hym2 = 1./hy/hy; hzm2 = 1./hz/hz; 1572 | 1573 | /* Get local grid boundaries */ 1574 | ierr = DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm); CHKERRQ (ierr); 1575 | ierr = DAGetGhostCorners(user->da,&gxs,&gys,&gzs,&gxm,&gym,&gzm); 1576 | CHKERRQ (ierr); 1577 | 1578 | /* Determine the interior of the locally owned part of the grid. */ 1579 | if (user->period == DA_XYZPERIODIC) { 1580 | xints = xs; xinte = xs+xm; 1581 | yints = ys; yinte = ys+ym; 1582 | zints = zs; zinte = zs+zm; } 1583 | else { 1584 | SETERRQ(PETSC_ERR_ARG_WRONGSTATE, 1585 | "Only periodic boundary conditions in 3-D Cahn-Hilliard"); } 1586 | 1587 | /* Get parameters for matrix creation */ 1588 | i = xm * ym * zm * user->mc; 1589 | j = mx * my * mz * user->mc; 1590 | 1591 | /* Create the distributed alpha matrix */ 1592 | ierr = MatCreateMPIAIJ (user->comm, i,i, j,j, 25,PETSC_NULL, 12,PETSC_NULL, 1593 | &(user->alpha)); 1594 | CHKERRQ (ierr); 1595 | 1596 | /* Get the local-to-global mapping and associate it with the matrix */ 1597 | ierr = DAGetISLocalToGlobalMapping (user->da, &user->isltog); CHKERRQ (ierr); 1598 | ierr = MatSetLocalToGlobalMapping (user->alpha, user->isltog); CHKERRQ (ierr); 1599 | 1600 | /* Pre-compute alpha term values (see ch_residual_2d above) 1601 | This should be the only place they're actually computed; they'll be 1602 | used below. */ 1603 | avalue[0] = avalue[24] = -alpha*hzm2*hzm2; 1604 | avalue[1] = avalue[5] = avalue[19] = avalue[23] = -2.*alpha*hym2*hzm2; 1605 | avalue[2] = avalue[4] = avalue[20] = avalue[22] = -2.*alpha*hxm2*hzm2; 1606 | avalue[3] = avalue[21] = 4.*alpha*hzm2*(hxm2+hym2+hzm2); 1607 | avalue[6] = avalue[18] = -alpha*hym2*hym2; 1608 | avalue[7] = avalue[9] = avalue[15] = avalue[17] = -2.*alpha*hxm2*hym2; 1609 | avalue[8] = avalue[16] = 4.*alpha*hym2*(hxm2+hym2+hzm2); 1610 | avalue[10] = avalue[14] = -alpha*hxm2*hxm2; 1611 | avalue[11] = avalue[13] = 4.*alpha*hxm2*(hxm2+hym2+hzm2); 1612 | avalue[12] = -alpha*(6.*hxm2*hxm2 + 6.*hym2*hym2 + 6.*hzm2*hzm2 1613 | + 8.*hxm2*hym2 + 8.*hxm2*hzm2 + 8.*hym2*hzm2); 1614 | 1615 | /* Loop over interior nodes */ 1616 | for (k=zints-gzs; k<zinte-gzs; k++) 1617 | for (j=k*gym + yints-gys; j<k*gym + yinte-gys; j++) 1618 | for (i=j*gxm + xints-gxs; i<j*gxm + xinte-gxs; i++) { 1619 | 1620 | /* Form the column indices */ 1621 | columns[0]=C(i+2*gxm*gym); 1622 | columns[1]=C(i+gxm*gym+gxm); 1623 | columns[2]=C(i+gxm*gym-1); 1624 | columns[3]=C(i+gxm*gym); 1625 | columns[4]=C(i+gxm*gym+1); 1626 | columns[5]=C(i+gxm*gym-gxm); 1627 | columns[6]=C(i+2*gxm); 1628 | columns[7]=C(i+gxm-1); columns[8]=C(i+gxm); columns[9]=C(i+gxm+1); 1629 | columns[10]=C(i-2); columns[11]=C(i-1); columns[12]=C(i); 1630 | columns[13]=C(i+1); columns[14]=C(i+2); 1631 | columns[15]=C(i-gxm-1); columns[16]=C(i-gxm); columns[17]=C(i-gxm+1); 1632 | columns[18]=C(i-2*gxm); 1633 | columns[19]=C(i-gxm*gym+gxm); 1634 | columns[20]=C(i-gxm*gym-1); 1635 | columns[21]=C(i-gxm*gym); 1636 | columns[22]=C(i-gxm*gym+1); 1637 | columns[23]=C(i-gxm*gym-gxm); 1638 | columns[24]=C(i-2*gxm*gym); 1639 | 1640 | node = C(i); 1641 | MatSetValuesLocal (user->alpha, 1,&node, 25,columns, avalue, 1642 | INSERT_VALUES); 1643 | } 1644 | 1645 | /* Assemble matrix, using the 2-step process: 1646 | MatAssemblyBegin(), MatAssemblyEnd(). */ 1647 | ierr = MatAssemblyBegin(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 1648 | ierr = MatAssemblyEnd(user->alpha,MAT_FINAL_ASSEMBLY); CHKERRQ (ierr); 1649 | 1650 | /* Make J a copy of alpha with the same local mapping (setting new mapping is 1651 | unnecessary with PETSc 2.1.5). */ 1652 | ierr = MatDuplicate (user->alpha, MAT_COPY_VALUES, &(user->J)); CHKERRQ (ierr); 1653 | 1654 | /* Tell the matrix we will never add a new nonzero location to the 1655 | matrix. If we do, it will generate an error. */ 1656 | ierr = MatSetOption(user->J,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ (ierr); 1657 | 1658 | return ierr; 1659 | }