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