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