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