1! Program usage: mpiexec -n <proc> plate2f [all TAO options] 2! 3! This example demonstrates use of the TAO package to solve a bound constrained 4! minimization problem. This example is based on a problem from the 5! MINPACK-2 test suite. Given a rectangular 2-D domain and boundary values 6! along the edges of the domain, the objective is to find the surface 7! with the minimal area that satisfies the boundary conditions. 8! The command line options are: 9! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction 10! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction 11! -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction 12! -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction 13! -bheight <ht>, where <ht> = height of the plate 14! 15#include "petsc/finclude/petscdmda.h" 16#include "petsc/finclude/petsctao.h" 17module plate2fmodule 18 use petscdmda 19 use petsctao 20 21 implicit none 22 Vec localX, localV 23 Vec Top, Left 24 Vec Right, Bottom 25 DM dm 26 PetscReal bheight 27 PetscInt bmx, bmy 28! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 29! Variable declarations 30! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 31! 32! Variables: 33! (common from plate2f.h): 34! Nx, Ny number of processors in x- and y- directions 35! mx, my number of grid points in x,y directions 36! N global dimension of vector 37 PetscInt mx, my, Nx, Ny, N 38contains 39 40! --------------------------------------------------------------------- 41! 42! FormFunctionGradient - Evaluates function f(X). 43! 44! Input Parameters: 45! tao - the Tao context 46! X - the input vector 47! dummy - optional user-defined context, as set by TaoSetFunction() 48! (not used here) 49! 50! Output Parameters: 51! fcn - the newly evaluated function 52! G - the gradient vector 53! info - error code 54! 55 56 subroutine FormFunctionGradient(ta, X, fcn, G, dummy, ierr) 57! Input/output variables 58 59 Tao ta 60 PetscReal fcn 61 Vec X, G 62 PetscErrorCode ierr 63 PetscInt dummy 64 65 PetscInt i, j, row 66 PetscInt xs, xm 67 PetscInt gxs, gxm 68 PetscInt ys, ym 69 PetscInt gys, gym 70 PetscReal ft, zero, hx, hy, hydhx, hxdhy 71 PetscReal area, rhx, rhy 72 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3 73 PetscReal d4, d5, d6, d7, d8 74 PetscReal df1dxc, df2dxc, df3dxc, df4dxc 75 PetscReal df5dxc, df6dxc 76 PetscReal xc, xl, xr, xt, xb, xlt, xrb 77 78 PetscReal, pointer :: g_v(:), x_v(:) 79 PetscReal, pointer :: top_v(:), left_v(:) 80 PetscReal, pointer :: right_v(:), bottom_v(:) 81 82 ft = 0.0 83 zero = 0.0 84 hx = 1.0/real(mx + 1) 85 hy = 1.0/real(my + 1) 86 hydhx = hy/hx 87 hxdhy = hx/hy 88 area = 0.5*hx*hy 89 rhx = real(mx) + 1.0 90 rhy = real(my) + 1.0 91 92! Get local mesh boundaries 93 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 94 PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)) 95 96! Scatter ghost points to local vector 97 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr)) 98 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr)) 99 100! Initialize the vector to zero 101 PetscCall(VecSet(localV, zero, ierr)) 102 103! Get arrays to vector data (See note above about using VecGetArray in Fortran) 104 PetscCall(VecGetArray(localX, x_v, ierr)) 105 PetscCall(VecGetArray(localV, g_v, ierr)) 106 PetscCall(VecGetArray(Top, top_v, ierr)) 107 PetscCall(VecGetArray(Bottom, bottom_v, ierr)) 108 PetscCall(VecGetArray(Left, left_v, ierr)) 109 PetscCall(VecGetArray(Right, right_v, ierr)) 110 111! Compute function over the locally owned part of the mesh 112 do j = ys, ys + ym - 1 113 do i = xs, xs + xm - 1 114 row = (j - gys)*gxm + (i - gxs) 115 xc = x_v(1 + row) 116 xt = xc 117 xb = xc 118 xr = xc 119 xl = xc 120 xrb = xc 121 xlt = xc 122 123 if (i == 0) then !left side 124 xl = left_v(1 + j - ys + 1) 125 xlt = left_v(1 + j - ys + 2) 126 else 127 xl = x_v(1 + row - 1) 128 end if 129 130 if (j == 0) then !bottom side 131 xb = bottom_v(1 + i - xs + 1) 132 xrb = bottom_v(1 + i - xs + 2) 133 else 134 xb = x_v(1 + row - gxm) 135 end if 136 137 if (i + 1 == gxs + gxm) then !right side 138 xr = right_v(1 + j - ys + 1) 139 xrb = right_v(1 + j - ys) 140 else 141 xr = x_v(1 + row + 1) 142 end if 143 144 if (j + 1 == gys + gym) then !top side 145 xt = top_v(1 + i - xs + 1) 146 xlt = top_v(1 + i - xs) 147 else 148 xt = x_v(1 + row + gxm) 149 end if 150 151 if ((i > gxs) .and. (j + 1 < gys + gym)) then 152 xlt = x_v(1 + row - 1 + gxm) 153 end if 154 155 if ((j > gys) .and. (i + 1 < gxs + gxm)) then 156 xrb = x_v(1 + row + 1 - gxm) 157 end if 158 159 d1 = xc - xl 160 d2 = xc - xr 161 d3 = xc - xt 162 d4 = xc - xb 163 d5 = xr - xrb 164 d6 = xrb - xb 165 d7 = xlt - xl 166 d8 = xt - xlt 167 168 df1dxc = d1*hydhx 169 df2dxc = d1*hydhx + d4*hxdhy 170 df3dxc = d3*hxdhy 171 df4dxc = d2*hydhx + d3*hxdhy 172 df5dxc = d2*hydhx 173 df6dxc = d4*hxdhy 174 175 d1 = d1*rhx 176 d2 = d2*rhx 177 d3 = d3*rhy 178 d4 = d4*rhy 179 d5 = d5*rhy 180 d6 = d6*rhx 181 d7 = d7*rhy 182 d8 = d8*rhx 183 184 f1 = sqrt(1.0 + d1*d1 + d7*d7) 185 f2 = sqrt(1.0 + d1*d1 + d4*d4) 186 f3 = sqrt(1.0 + d3*d3 + d8*d8) 187 f4 = sqrt(1.0 + d3*d3 + d2*d2) 188 f5 = sqrt(1.0 + d2*d2 + d5*d5) 189 f6 = sqrt(1.0 + d4*d4 + d6*d6) 190 191 ft = ft + f2 + f4 192 193 df1dxc = df1dxc/f1 194 df2dxc = df2dxc/f2 195 df3dxc = df3dxc/f3 196 df4dxc = df4dxc/f4 197 df5dxc = df5dxc/f5 198 df6dxc = df6dxc/f6 199 200 g_v(1 + row) = 0.5*(df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) 201 end do 202 end do 203 204! Compute triangular areas along the border of the domain. 205 if (xs == 0) then ! left side 206 do j = ys, ys + ym - 1 207 d3 = (left_v(1 + j - ys + 1) - left_v(1 + j - ys + 2))*rhy 208 d2 = (left_v(1 + j - ys + 1) - x_v(1 + (j - gys)*gxm))*rhx 209 ft = ft + sqrt(1.0 + d3*d3 + d2*d2) 210 end do 211 end if 212 213 if (ys == 0) then !bottom side 214 do i = xs, xs + xm - 1 215 d2 = (bottom_v(1 + i + 1 - xs) - bottom_v(1 + i - xs + 2))*rhx 216 d3 = (bottom_v(1 + i - xs + 1) - x_v(1 + i - gxs))*rhy 217 ft = ft + sqrt(1.0 + d3*d3 + d2*d2) 218 end do 219 end if 220 221 if (xs + xm == mx) then ! right side 222 do j = ys, ys + ym - 1 223 d1 = (x_v(1 + (j + 1 - gys)*gxm - 1) - right_v(1 + j - ys + 1))*rhx 224 d4 = (right_v(1 + j - ys) - right_v(1 + j - ys + 1))*rhy 225 ft = ft + sqrt(1.0 + d1*d1 + d4*d4) 226 end do 227 end if 228 229 if (ys + ym == my) then 230 do i = xs, xs + xm - 1 231 d1 = (x_v(1 + (gym - 1)*gxm + i - gxs) - top_v(1 + i - xs + 1))*rhy 232 d4 = (top_v(1 + i - xs + 1) - top_v(1 + i - xs))*rhx 233 ft = ft + sqrt(1.0 + d1*d1 + d4*d4) 234 end do 235 end if 236 237 if ((ys == 0) .and. (xs == 0)) then 238 d1 = (left_v(1 + 0) - left_v(1 + 1))*rhy 239 d2 = (bottom_v(1 + 0) - bottom_v(1 + 1))*rhx 240 ft = ft + sqrt(1.0 + d1*d1 + d2*d2) 241 end if 242 243 if ((ys + ym == my) .and. (xs + xm == mx)) then 244 d1 = (right_v(1 + ym + 1) - right_v(1 + ym))*rhy 245 d2 = (top_v(1 + xm + 1) - top_v(1 + xm))*rhx 246 ft = ft + sqrt(1.0 + d1*d1 + d2*d2) 247 end if 248 249 ft = ft*area 250 PetscCallMPI(MPI_Allreduce(ft, fcn, 1, MPIU_REAL, MPIU_SUM, PETSC_COMM_WORLD, ierr)) 251 252! Restore vectors 253 PetscCall(VecRestoreArray(localX, x_v, ierr)) 254 PetscCall(VecRestoreArray(localV, g_v, ierr)) 255 PetscCall(VecRestoreArray(Left, left_v, ierr)) 256 PetscCall(VecRestoreArray(Top, top_v, ierr)) 257 PetscCall(VecRestoreArray(Bottom, bottom_v, ierr)) 258 PetscCall(VecRestoreArray(Right, right_v, ierr)) 259 260! Scatter values to global vector 261 PetscCall(DMLocalToGlobalBegin(dm, localV, INSERT_VALUES, G, ierr)) 262 PetscCall(DMLocalToGlobalEnd(dm, localV, INSERT_VALUES, G, ierr)) 263 264 PetscCall(PetscLogFlops(70.0d0*xm*ym, ierr)) 265 266 end !FormFunctionGradient 267 268! ---------------------------------------------------------------------------- 269! 270! 271! FormHessian - Evaluates Hessian matrix. 272! 273! Input Parameters: 274!. tao - the Tao context 275!. X - input vector 276!. dummy - not used 277! 278! Output Parameters: 279!. Hessian - Hessian matrix 280!. Hpc - optionally different matrix used to construct the preconditioner 281! 282! Notes: 283! Due to mesh point reordering with DMs, we must always work 284! with the local mesh points, and then transform them to the new 285! global numbering with the local-to-global mapping. We cannot work 286! directly with the global numbers for the original uniprocessor mesh! 287! 288! MatSetValuesLocal(), using the local ordering (including 289! ghost points!) 290! - Then set matrix entries using the local ordering 291! by calling MatSetValuesLocal() 292 293 subroutine FormHessian(ta, X, Hessian, Hpc, dummy, ierr) 294 295 Tao ta 296 Vec X 297 Mat Hessian, Hpc 298 PetscInt dummy 299 PetscErrorCode ierr 300 301 PetscInt i, j, k, row 302 PetscInt xs, xm, gxs, gxm 303 PetscInt ys, ym, gys, gym 304 PetscInt col(0:6) 305 PetscReal hx, hy, hydhx, hxdhy, rhx, rhy 306 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3 307 PetscReal d4, d5, d6, d7, d8 308 PetscReal xc, xl, xr, xt, xb, xlt, xrb 309 PetscReal hl, hr, ht, hb, hc, htl, hbr 310 311 PetscReal, pointer :: right_v(:), left_v(:) 312 PetscReal, pointer :: bottom_v(:), top_v(:) 313 PetscReal, pointer :: x_v(:) 314 PetscReal v(0:6) 315 PetscBool assembled 316 PetscInt i1 317 318 i1 = 1 319 320! Set various matrix options 321 PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE, ierr)) 322 323! Get local mesh boundaries 324 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 325 PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)) 326 327! Scatter ghost points to local vectors 328 PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX, ierr)) 329 PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX, ierr)) 330 331! Get pointers to vector data (see note on Fortran arrays above) 332 PetscCall(VecGetArray(localX, x_v, ierr)) 333 PetscCall(VecGetArray(Top, top_v, ierr)) 334 PetscCall(VecGetArray(Bottom, bottom_v, ierr)) 335 PetscCall(VecGetArray(Left, left_v, ierr)) 336 PetscCall(VecGetArray(Right, right_v, ierr)) 337 338! Initialize matrix entries to zero 339 PetscCall(MatAssembled(Hessian, assembled, ierr)) 340 if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian, ierr)) 341 342 rhx = real(mx) + 1.0 343 rhy = real(my) + 1.0 344 hx = 1.0/rhx 345 hy = 1.0/rhy 346 hydhx = hy/hx 347 hxdhy = hx/hy 348! compute Hessian over the locally owned part of the mesh 349 350 do i = xs, xs + xm - 1 351 do j = ys, ys + ym - 1 352 row = (j - gys)*gxm + (i - gxs) 353 354 xc = x_v(1 + row) 355 xt = xc 356 xb = xc 357 xr = xc 358 xl = xc 359 xrb = xc 360 xlt = xc 361 362 if (i == gxs) then ! Left side 363 xl = left_v(1 + j - ys + 1) 364 xlt = left_v(1 + j - ys + 2) 365 else 366 xl = x_v(1 + row - 1) 367 end if 368 369 if (j == gys) then ! bottom side 370 xb = bottom_v(1 + i - xs + 1) 371 xrb = bottom_v(1 + i - xs + 2) 372 else 373 xb = x_v(1 + row - gxm) 374 end if 375 376 if (i + 1 == gxs + gxm) then !right side 377 xr = right_v(1 + j - ys + 1) 378 xrb = right_v(1 + j - ys) 379 else 380 xr = x_v(1 + row + 1) 381 end if 382 383 if (j + 1 == gym + gys) then !top side 384 xt = top_v(1 + i - xs + 1) 385 xlt = top_v(1 + i - xs) 386 else 387 xt = x_v(1 + row + gxm) 388 end if 389 390 if ((i > gxs) .and. (j + 1 < gys + gym)) then 391 xlt = x_v(1 + row - 1 + gxm) 392 end if 393 394 if ((i + 1 < gxs + gxm) .and. (j > gys)) then 395 xrb = x_v(1 + row + 1 - gxm) 396 end if 397 398 d1 = (xc - xl)*rhx 399 d2 = (xc - xr)*rhx 400 d3 = (xc - xt)*rhy 401 d4 = (xc - xb)*rhy 402 d5 = (xrb - xr)*rhy 403 d6 = (xrb - xb)*rhx 404 d7 = (xlt - xl)*rhy 405 d8 = (xlt - xt)*rhx 406 407 f1 = sqrt(1.0 + d1*d1 + d7*d7) 408 f2 = sqrt(1.0 + d1*d1 + d4*d4) 409 f3 = sqrt(1.0 + d3*d3 + d8*d8) 410 f4 = sqrt(1.0 + d3*d3 + d2*d2) 411 f5 = sqrt(1.0 + d2*d2 + d5*d5) 412 f6 = sqrt(1.0 + d4*d4 + d6*d6) 413 414 hl = (-hydhx*(1.0 + d7*d7) + d1*d7)/(f1*f1*f1) + (-hydhx*(1.0 + d4*d4) + d1*d4)/(f2*f2*f2) 415 416 hr = (-hydhx*(1.0 + d5*d5) + d2*d5)/(f5*f5*f5) + (-hydhx*(1.0 + d3*d3) + d2*d3)/(f4*f4*f4) 417 418 ht = (-hxdhy*(1.0 + d8*d8) + d3*d8)/(f3*f3*f3) + (-hxdhy*(1.0 + d2*d2) + d2*d3)/(f4*f4*f4) 419 420 hb = (-hxdhy*(1.0 + d6*d6) + d4*d6)/(f6*f6*f6) + (-hxdhy*(1.0 + d1*d1) + d1*d4)/(f2*f2*f2) 421 422 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6) 423 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3) 424 425 hc = hydhx*(1.0 + d7*d7)/(f1*f1*f1) + hxdhy*(1.0 + d8*d8)/(f3*f3*f3) + hydhx*(1.0 + d5*d5)/(f5*f5*f5) + & 426 & hxdhy*(1.0 + d6*d6)/(f6*f6*f6) + (hxdhy*(1.0 + d1*d1) + hydhx*(1.0 + d4*d4) - 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0 + d2*d2) + & 427 & hydhx*(1.0 + d3*d3) - 2*d2*d3)/(f4*f4*f4) 428 429 hl = hl*0.5 430 hr = hr*0.5 431 ht = ht*0.5 432 hb = hb*0.5 433 hbr = hbr*0.5 434 htl = htl*0.5 435 hc = hc*0.5 436 437 k = 0 438 439 if (j > 0) then 440 v(k) = hb 441 col(k) = row - gxm 442 k = k + 1 443 end if 444 445 if ((j > 0) .and. (i < mx - 1)) then 446 v(k) = hbr 447 col(k) = row - gxm + 1 448 k = k + 1 449 end if 450 451 if (i > 0) then 452 v(k) = hl 453 col(k) = row - 1 454 k = k + 1 455 end if 456 457 v(k) = hc 458 col(k) = row 459 k = k + 1 460 461 if (i < mx - 1) then 462 v(k) = hr 463 col(k) = row + 1 464 k = k + 1 465 end if 466 467 if ((i > 0) .and. (j < my - 1)) then 468 v(k) = htl 469 col(k) = row + gxm - 1 470 k = k + 1 471 end if 472 473 if (j < my - 1) then 474 v(k) = ht 475 col(k) = row + gxm 476 k = k + 1 477 end if 478 479! Set matrix values using local numbering, defined earlier in main routine 480 PetscCall(MatSetValuesLocal(Hessian, i1, [row], k, col, v, INSERT_VALUES, ierr)) 481 482 end do 483 end do 484 485! restore vectors 486 PetscCall(VecRestoreArray(localX, x_v, ierr)) 487 PetscCall(VecRestoreArray(Left, left_v, ierr)) 488 PetscCall(VecRestoreArray(Right, right_v, ierr)) 489 PetscCall(VecRestoreArray(Top, top_v, ierr)) 490 PetscCall(VecRestoreArray(Bottom, bottom_v, ierr)) 491 492! Assemble the matrix 493 PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY, ierr)) 494 PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY, ierr)) 495 496 PetscCall(PetscLogFlops(199.0d0*xm*ym, ierr)) 497 498 end 499 500! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h 501 502! ---------------------------------------------------------------------------- 503! 504!/* 505! MSA_BoundaryConditions - calculates the boundary conditions for the region 506! 507! 508!*/ 509 510 subroutine MSA_BoundaryConditions(ierr) 511 512 PetscErrorCode ierr 513 PetscInt i, j, k, limit, maxits 514 PetscInt xs, xm, gxs, gxm 515 PetscInt ys, ym, gys, gym 516 PetscInt bsize, lsize 517 PetscInt tsize, rsize 518 PetscReal one, two, three, tol 519 PetscReal scl, fnorm, det, xt 520 PetscReal yt, hx, hy, u1, u2, nf1, nf2 521 PetscReal njac11, njac12, njac21, njac22 522 PetscReal b, t, l, r 523 PetscReal, pointer :: boundary_v(:) 524 525 logical exitloop 526 PetscBool flg 527 528 limit = 0 529 maxits = 5 530 tol = 1e-10 531 b = -0.5 532 t = 0.5 533 l = -0.5 534 r = 0.5 535 xt = 0 536 yt = 0 537 one = 1.0 538 two = 2.0 539 three = 3.0 540 541 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 542 PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)) 543 544 bsize = xm + 2 545 lsize = ym + 2 546 rsize = ym + 2 547 tsize = xm + 2 548 549 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, bsize, PETSC_DECIDE, Bottom, ierr)) 550 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, tsize, PETSC_DECIDE, Top, ierr)) 551 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, lsize, PETSC_DECIDE, Left, ierr)) 552 PetscCall(VecCreateMPI(PETSC_COMM_WORLD, rsize, PETSC_DECIDE, Right, ierr)) 553 554 hx = (r - l)/(mx + 1) 555 hy = (t - b)/(my + 1) 556 557 do j = 0, 3 558 559 if (j == 0) then 560 yt = b 561 xt = l + hx*xs 562 limit = bsize 563 PetscCall(VecGetArray(Bottom, boundary_v, ierr)) 564 565 elseif (j == 1) then 566 yt = t 567 xt = l + hx*xs 568 limit = tsize 569 PetscCall(VecGetArray(Top, boundary_v, ierr)) 570 571 elseif (j == 2) then 572 yt = b + hy*ys 573 xt = l 574 limit = lsize 575 PetscCall(VecGetArray(Left, boundary_v, ierr)) 576 577 elseif (j == 3) then 578 yt = b + hy*ys 579 xt = r 580 limit = rsize 581 PetscCall(VecGetArray(Right, boundary_v, ierr)) 582 end if 583 584 do i = 0, limit - 1 585 586 u1 = xt 587 u2 = -yt 588 k = 0 589 exitloop = .false. 590 do while (k < maxits .and. (.not. exitloop)) 591 592 nf1 = u1 + u1*u2*u2 - u1*u1*u1/three - xt 593 nf2 = -u2 - u1*u1*u2 + u2*u2*u2/three - yt 594 fnorm = sqrt(nf1*nf1 + nf2*nf2) 595 if (fnorm > tol) then 596 njac11 = one + u2*u2 - u1*u1 597 njac12 = two*u1*u2 598 njac21 = -two*u1*u2 599 njac22 = -one - u1*u1 + u2*u2 600 det = njac11*njac22 - njac21*njac12 601 u1 = u1 - (njac22*nf1 - njac12*nf2)/det 602 u2 = u2 - (njac11*nf2 - njac21*nf1)/det 603 else 604 exitloop = .true. 605 end if 606 k = k + 1 607 end do 608 609 boundary_v(1 + i) = u1*u1 - u2*u2 610 if ((j == 0) .or. (j == 1)) then 611 xt = xt + hx 612 else 613 yt = yt + hy 614 end if 615 616 end do 617 618 if (j == 0) then 619 PetscCall(VecRestoreArray(Bottom, boundary_v, ierr)) 620 elseif (j == 1) then 621 PetscCall(VecRestoreArray(Top, boundary_v, ierr)) 622 elseif (j == 2) then 623 PetscCall(VecRestoreArray(Left, boundary_v, ierr)) 624 elseif (j == 3) then 625 PetscCall(VecRestoreArray(Right, boundary_v, ierr)) 626 end if 627 628 end do 629 630! Scale the boundary if desired 631 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bottom', scl, flg, ierr)) 632 if (flg .eqv. PETSC_TRUE) then 633 PetscCall(VecScale(Bottom, scl, ierr)) 634 end if 635 636 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-top', scl, flg, ierr)) 637 if (flg .eqv. PETSC_TRUE) then 638 PetscCall(VecScale(Top, scl, ierr)) 639 end if 640 641 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-right', scl, flg, ierr)) 642 if (flg .eqv. PETSC_TRUE) then 643 PetscCall(VecScale(Right, scl, ierr)) 644 end if 645 646 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-left', scl, flg, ierr)) 647 if (flg .eqv. PETSC_TRUE) then 648 PetscCall(VecScale(Left, scl, ierr)) 649 end if 650 651 end 652 653! ---------------------------------------------------------------------------- 654! 655!/* 656! MSA_Plate - Calculates an obstacle for surface to stretch over 657! 658! Output Parameter: 659!. xl - lower bound vector 660!. xu - upper bound vector 661! 662!*/ 663 664 subroutine MSA_Plate(ta, xl, xu, dummy, ierr) 665 666 Tao ta 667 Vec xl, xu 668 PetscErrorCode ierr 669 PetscInt i, j, row 670 PetscInt xs, xm, ys, ym 671 PetscReal lb, ub 672 PetscInt dummy 673 PetscReal, pointer :: xl_v(:) 674 675 lb = PETSC_NINFINITY 676 ub = PETSC_INFINITY 677 678 if (bmy < 0) bmy = 0 679 if (bmy > my) bmy = my 680 if (bmx < 0) bmx = 0 681 if (bmx > mx) bmx = mx 682 683 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 684 685 PetscCall(VecSet(xl, lb, ierr)) 686 PetscCall(VecSet(xu, ub, ierr)) 687 688 PetscCall(VecGetArray(xl, xl_v, ierr)) 689 690 do i = xs, xs + xm - 1 691 692 do j = ys, ys + ym - 1 693 694 row = (j - ys)*xm + (i - xs) 695 696 if (i >= ((mx - bmx)/2) .and. i < (mx - (mx - bmx)/2) .and. & 697 & j >= ((my - bmy)/2) .and. j < (my - (my - bmy)/2)) then 698 xl_v(1 + row) = bheight 699 700 end if 701 702 end do 703 end do 704 705 PetscCall(VecRestoreArray(xl, xl_v, ierr)) 706 707 end 708 709! ---------------------------------------------------------------------------- 710! 711!/* 712! MSA_InitialPoint - Calculates an obstacle for surface to stretch over 713! 714! Output Parameter: 715!. X - vector for initial guess 716! 717!*/ 718 719 subroutine MSA_InitialPoint(X, ierr) 720 721 Vec X 722 PetscErrorCode ierr 723 PetscInt start, i, j 724 PetscInt row 725 PetscInt xs, xm, gxs, gxm 726 PetscInt ys, ym, gys, gym 727 PetscReal zero, np5 728 729 PetscReal, pointer :: left_v(:), right_v(:) 730 PetscReal, pointer :: bottom_v(:), top_v(:) 731 PetscReal, pointer :: x_v(:) 732 PetscBool flg 733 PetscRandom rctx 734 735 zero = 0.0 736 np5 = -0.5 737 738 PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-start', start, flg, ierr)) 739 740 if ((flg .eqv. PETSC_TRUE) .and. (start == 0)) then ! the zero vector is reasonable 741 PetscCall(VecSet(X, zero, ierr)) 742 743 elseif ((flg .eqv. PETSC_TRUE) .and. (start > 0)) then ! random start -0.5 < xi < 0.5 744 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, rctx, ierr)) 745 do i = 0, start - 1 746 PetscCall(VecSetRandom(X, rctx, ierr)) 747 end do 748 749 PetscCall(PetscRandomDestroy(rctx, ierr)) 750 PetscCall(VecShift(X, np5, ierr)) 751 752 else ! average of boundary conditions 753 754! Get Local mesh boundaries 755 PetscCall(DMDAGetCorners(dm, xs, ys, PETSC_NULL_INTEGER, xm, ym, PETSC_NULL_INTEGER, ierr)) 756 PetscCall(DMDAGetGhostCorners(dm, gxs, gys, PETSC_NULL_INTEGER, gxm, gym, PETSC_NULL_INTEGER, ierr)) 757 758! Get pointers to vector data 759 PetscCall(VecGetArray(Top, top_v, ierr)) 760 PetscCall(VecGetArray(Bottom, bottom_v, ierr)) 761 PetscCall(VecGetArray(Left, left_v, ierr)) 762 PetscCall(VecGetArray(Right, right_v, ierr)) 763 764 PetscCall(VecGetArray(localX, x_v, ierr)) 765 766! Perform local computations 767 do j = ys, ys + ym - 1 768 do i = xs, xs + xm - 1 769 row = (j - gys)*gxm + (i - gxs) 770 x_v(1 + row) = ((j + 1)*bottom_v(1 + i - xs + 1)/my + (my - j + 1)*top_v(1 + i - xs + 1)/(my + 2) + & 771 & (i + 1)*left_v(1 + j - ys + 1)/mx + (mx - i + 1)*right_v(1 + j - ys + 1)/(mx + 2))*0.5 772 end do 773 end do 774 775! Restore vectors 776 PetscCall(VecRestoreArray(localX, x_v, ierr)) 777 778 PetscCall(VecRestoreArray(Left, left_v, ierr)) 779 PetscCall(VecRestoreArray(Top, top_v, ierr)) 780 PetscCall(VecRestoreArray(Bottom, bottom_v, ierr)) 781 PetscCall(VecRestoreArray(Right, right_v, ierr)) 782 783 PetscCall(DMLocalToGlobalBegin(dm, localX, INSERT_VALUES, X, ierr)) 784 PetscCall(DMLocalToGlobalEnd(dm, localX, INSERT_VALUES, X, ierr)) 785 786 end if 787 788 end 789 790end module 791 792program plate2f 793 use plate2fmodule 794 implicit none 795 796 PetscErrorCode ierr ! used to check for functions returning nonzeros 797 Vec x ! solution vector 798 PetscInt m ! number of local elements in vector 799 Tao ta ! Tao solver context 800 Mat H ! Hessian matrix 801 ISLocalToGlobalMapping isltog ! local to global mapping object 802 PetscBool flg 803 PetscInt i1, i3, i7 804 805! Initialize Tao 806 807 i1 = 1 808 i3 = 3 809 i7 = 7 810 811 PetscCallA(PetscInitialize(ierr)) 812 813! Specify default dimensions of the problem 814 mx = 10 815 my = 10 816 bheight = 0.1 817 818! Check for any command line arguments that override defaults 819 820 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr)) 821 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr)) 822 823 bmx = mx/2 824 bmy = my/2 825 826 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmx', bmx, flg, ierr)) 827 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bmy', bmy, flg, ierr)) 828 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-bheight', bheight, flg, ierr)) 829 830! Calculate any derived values from parameters 831 N = mx*my 832 833! Let PETSc determine the dimensions of the local vectors 834 Nx = PETSC_DECIDE 835 NY = PETSC_DECIDE 836 837! A two dimensional distributed array will help define this problem, which 838! derives from an elliptic PDE on a two-dimensional domain. From the 839! distributed array, create the vectors 840 841 PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, Nx, Ny, i1, i1, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm, ierr)) 842 PetscCallA(DMSetFromOptions(dm, ierr)) 843 PetscCallA(DMSetUp(dm, ierr)) 844 845! Extract global and local vectors from DM; The local vectors are 846! used solely as work space for the evaluation of the function, 847! gradient, and Hessian. Duplicate for remaining vectors that are 848! the same types. 849 850 PetscCallA(DMCreateGlobalVector(dm, x, ierr)) 851 PetscCallA(DMCreateLocalVector(dm, localX, ierr)) 852 PetscCallA(VecDuplicate(localX, localV, ierr)) 853 854! Create a matrix data structure to store the Hessian. 855! Here we (optionally) also associate the local numbering scheme 856! with the matrix so that later we can use local indices for matrix 857! assembly 858 859 PetscCallA(VecGetLocalSize(x, m, ierr)) 860 PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, m, m, N, N, i7, PETSC_NULL_INTEGER_ARRAY, i3, PETSC_NULL_INTEGER_ARRAY, H, ierr)) 861 862 PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr)) 863 PetscCallA(DMGetLocalToGlobalMapping(dm, isltog, ierr)) 864 PetscCallA(MatSetLocalToGlobalMapping(H, isltog, isltog, ierr)) 865 866! The Tao code begins here 867! Create TAO solver and set desired solution method. 868! This problems uses bounded variables, so the 869! method must either be 'tao_tron' or 'tao_blmvm' 870 871 PetscCallA(TaoCreate(PETSC_COMM_WORLD, ta, ierr)) 872 PetscCallA(TaoSetType(ta, TAOBLMVM, ierr)) 873 874! Set minimization function and gradient, hessian evaluation functions 875 876 PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr)) 877 878 PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr)) 879 880! Set Variable bounds 881 PetscCallA(MSA_BoundaryConditions(ierr)) 882 PetscCallA(TaoSetVariableBoundsRoutine(ta, MSA_Plate, 0, ierr)) 883 884! Set the initial solution guess 885 PetscCallA(MSA_InitialPoint(x, ierr)) 886 PetscCallA(TaoSetSolution(ta, x, ierr)) 887 888! Check for any tao command line options 889 PetscCallA(TaoSetFromOptions(ta, ierr)) 890 891! Solve the application 892 PetscCallA(TaoSolve(ta, ierr)) 893 894! Free TAO data structures 895 PetscCallA(TaoDestroy(ta, ierr)) 896 897! Free PETSc data structures 898 PetscCallA(VecDestroy(x, ierr)) 899 PetscCallA(VecDestroy(Top, ierr)) 900 PetscCallA(VecDestroy(Bottom, ierr)) 901 PetscCallA(VecDestroy(Left, ierr)) 902 PetscCallA(VecDestroy(Right, ierr)) 903 PetscCallA(MatDestroy(H, ierr)) 904 PetscCallA(VecDestroy(localX, ierr)) 905 PetscCallA(VecDestroy(localV, ierr)) 906 PetscCallA(DMDestroy(dm, ierr)) 907 908! Finalize TAO 909 910 PetscCallA(PetscFinalize(ierr)) 911 912end program plate2f 913! 914!/*TEST 915! 916! build: 917! requires: !complex 918! 919! test: 920! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4 921! filter: sort -b 922! filter_output: sort -b 923! requires: !single 924! 925! test: 926! suffix: 2 927! nsize: 2 928! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4 929! filter: sort -b 930! filter_output: sort -b 931! requires: !single 932! 933!TEST*/ 934