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