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