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 PetscReal, pointer :: g_v(:),x_v(:) 208 PetscReal, pointer :: top_v(:),left_v(:) 209 PetscReal, pointer :: right_v(:),bottom_v(:) 210 211 ft = 0.0 212 zero = 0.0 213 hx = 1.0/real(mx + 1) 214 hy = 1.0/real(my + 1) 215 hydhx = hy/hx 216 hxdhy = hx/hy 217 area = 0.5 * hx * hy 218 rhx = real(mx) + 1.0 219 rhy = real(my) + 1.0 220 221! Get local mesh boundaries 222 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 223 PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)) 224 225! Scatter ghost points to local vector 226 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)) 227 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)) 228 229! Initialize the vector to zero 230 PetscCall(VecSet(localV,zero,ierr)) 231 232! Get arrays to vector data (See note above about using VecGetArrayF90 in Fortran) 233 PetscCall(VecGetArrayF90(localX,x_v,ierr)) 234 PetscCall(VecGetArrayF90(localV,g_v,ierr)) 235 PetscCall(VecGetArrayF90(Top,top_v,ierr)) 236 PetscCall(VecGetArrayF90(Bottom,bottom_v,ierr)) 237 PetscCall(VecGetArrayF90(Left,left_v,ierr)) 238 PetscCall(VecGetArrayF90(Right,right_v,ierr)) 239 240! Compute function over the locally owned part of the mesh 241 do j = ys,ys+ym-1 242 do i = xs,xs+xm-1 243 row = (j-gys)*gxm + (i-gxs) 244 xc = x_v(1+row) 245 xt = xc 246 xb = xc 247 xr = xc 248 xl = xc 249 xrb = xc 250 xlt = xc 251 252 if (i .eq. 0) then !left side 253 xl = left_v(1+j - ys + 1) 254 xlt = left_v(1+j - ys + 2) 255 else 256 xl = x_v(1+row - 1) 257 endif 258 259 if (j .eq. 0) then !bottom side 260 xb = bottom_v(1+i - xs + 1) 261 xrb = bottom_v(1+i - xs + 2) 262 else 263 xb = x_v(1+row - gxm) 264 endif 265 266 if (i + 1 .eq. gxs + gxm) then !right side 267 xr = right_v(1+j - ys + 1) 268 xrb = right_v(1+j - ys) 269 else 270 xr = x_v(1+row + 1) 271 endif 272 273 if (j + 1 .eq. gys + gym) then !top side 274 xt = top_v(1+i - xs + 1) 275 xlt = top_v(1+i - xs) 276 else 277 xt = x_v(1+row + gxm) 278 endif 279 280 if ((i .gt. gxs) .and. (j + 1 .lt. gys + gym)) then 281 xlt = x_v(1+row - 1 + gxm) 282 endif 283 284 if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then 285 xrb = x_v(1+row + 1 - gxm) 286 endif 287 288 d1 = xc-xl 289 d2 = xc-xr 290 d3 = xc-xt 291 d4 = xc-xb 292 d5 = xr-xrb 293 d6 = xrb-xb 294 d7 = xlt-xl 295 d8 = xt-xlt 296 297 df1dxc = d1 * hydhx 298 df2dxc = d1 * hydhx + d4 * hxdhy 299 df3dxc = d3 * hxdhy 300 df4dxc = d2 * hydhx + d3 * hxdhy 301 df5dxc = d2 * hydhx 302 df6dxc = d4 * hxdhy 303 304 d1 = d1 * rhx 305 d2 = d2 * rhx 306 d3 = d3 * rhy 307 d4 = d4 * rhy 308 d5 = d5 * rhy 309 d6 = d6 * rhx 310 d7 = d7 * rhy 311 d8 = d8 * rhx 312 313 f1 = sqrt(1.0 + d1*d1 + d7*d7) 314 f2 = sqrt(1.0 + d1*d1 + d4*d4) 315 f3 = sqrt(1.0 + d3*d3 + d8*d8) 316 f4 = sqrt(1.0 + d3*d3 + d2*d2) 317 f5 = sqrt(1.0 + d2*d2 + d5*d5) 318 f6 = sqrt(1.0 + d4*d4 + d6*d6) 319 320 ft = ft + f2 + f4 321 322 df1dxc = df1dxc / f1 323 df2dxc = df2dxc / f2 324 df3dxc = df3dxc / f3 325 df4dxc = df4dxc / f4 326 df5dxc = df5dxc / f5 327 df6dxc = df6dxc / f6 328 329 g_v(1+row) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) 330 enddo 331 enddo 332 333! Compute triangular areas along the border of the domain. 334 if (xs .eq. 0) then ! left side 335 do j=ys,ys+ym-1 336 d3 = (left_v(1+j-ys+1) - left_v(1+j-ys+2)) * rhy 337 d2 = (left_v(1+j-ys+1) - x_v(1+(j-gys)*gxm)) * rhx 338 ft = ft + sqrt(1.0 + d3*d3 + d2*d2) 339 enddo 340 endif 341 342 if (ys .eq. 0) then !bottom side 343 do i=xs,xs+xm-1 344 d2 = (bottom_v(1+i+1-xs)-bottom_v(1+i-xs+2)) * rhx 345 d3 = (bottom_v(1+i-xs+1)-x_v(1+i-gxs))*rhy 346 ft = ft + sqrt(1.0 + d3*d3 + d2*d2) 347 enddo 348 endif 349 350 if (xs + xm .eq. mx) then ! right side 351 do j=ys,ys+ym-1 352 d1 = (x_v(1+(j+1-gys)*gxm-1)-right_v(1+j-ys+1))*rhx 353 d4 = (right_v(1+j-ys) - right_v(1+j-ys+1))*rhy 354 ft = ft + sqrt(1.0 + d1*d1 + d4*d4) 355 enddo 356 endif 357 358 if (ys + ym .eq. my) then 359 do i=xs,xs+xm-1 360 d1 = (x_v(1+(gym-1)*gxm+i-gxs) - top_v(1+i-xs+1))*rhy 361 d4 = (top_v(1+i-xs+1) - top_v(1+i-xs))*rhx 362 ft = ft + sqrt(1.0 + d1*d1 + d4*d4) 363 enddo 364 endif 365 366 if ((ys .eq. 0) .and. (xs .eq. 0)) then 367 d1 = (left_v(1+0) - left_v(1+1)) * rhy 368 d2 = (bottom_v(1+0)-bottom_v(1+1))*rhx 369 ft = ft + sqrt(1.0 + d1*d1 + d2*d2) 370 endif 371 372 if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then 373 d1 = (right_v(1+ym+1) - right_v(1+ym))*rhy 374 d2 = (top_v(1+xm+1) - top_v(1+xm))*rhx 375 ft = ft + sqrt(1.0 + d1*d1 + d2*d2) 376 endif 377 378 ft = ft * area 379 PetscCallMPI(MPI_Allreduce(ft,fcn,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD,ierr)) 380 381! Restore vectors 382 PetscCall(VecRestoreArrayF90(localX,x_v,ierr)) 383 PetscCall(VecRestoreArrayF90(localV,g_v,ierr)) 384 PetscCall(VecRestoreArrayF90(Left,left_v,ierr)) 385 PetscCall(VecRestoreArrayF90(Top,top_v,ierr)) 386 PetscCall(VecRestoreArrayF90(Bottom,bottom_v,ierr)) 387 PetscCall(VecRestoreArrayF90(Right,right_v,ierr)) 388 389! Scatter values to global vector 390 PetscCall(DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr)) 391 PetscCall(DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr)) 392 393 PetscCall(PetscLogFlops(70.0d0*xm*ym,ierr)) 394 395 end !FormFunctionGradient 396 397! ---------------------------------------------------------------------------- 398! 399! 400! FormHessian - Evaluates Hessian matrix. 401! 402! Input Parameters: 403!. tao - the Tao context 404!. X - input vector 405!. dummy - not used 406! 407! Output Parameters: 408!. Hessian - Hessian matrix 409!. Hpc - optionally different preconditioning matrix 410!. flag - flag indicating matrix structure 411! 412! Notes: 413! Due to mesh point reordering with DMs, we must always work 414! with the local mesh points, and then transform them to the new 415! global numbering with the local-to-global mapping. We cannot work 416! directly with the global numbers for the original uniprocessor mesh! 417! 418! MatSetValuesLocal(), using the local ordering (including 419! ghost points!) 420! - Then set matrix entries using the local ordering 421! by calling MatSetValuesLocal() 422 423 subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr) 424 use plate2fmodule 425 implicit none 426 427 Tao tao 428 Vec X 429 Mat Hessian,Hpc 430 PetscInt dummy 431 PetscErrorCode ierr 432 433 PetscInt i,j,k,row 434 PetscInt xs,xm,gxs,gxm 435 PetscInt ys,ym,gys,gym 436 PetscInt col(0:6) 437 PetscReal hx,hy,hydhx,hxdhy,rhx,rhy 438 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3 439 PetscReal d4,d5,d6,d7,d8 440 PetscReal xc,xl,xr,xt,xb,xlt,xrb 441 PetscReal hl,hr,ht,hb,hc,htl,hbr 442 443 PetscReal,pointer :: right_v(:),left_v(:) 444 PetscReal,pointer :: bottom_v(:),top_v(:) 445 PetscReal,pointer :: x_v(:) 446 PetscReal v(0:6) 447 PetscBool assembled 448 PetscInt i1 449 450 i1=1 451 452! Set various matrix options 453 PetscCall(MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE,ierr)) 454 455! Get local mesh boundaries 456 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 457 PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)) 458 459! Scatter ghost points to local vectors 460 PetscCall(DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr)) 461 PetscCall(DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr)) 462 463! Get pointers to vector data (see note on Fortran arrays above) 464 PetscCall(VecGetArrayF90(localX,x_v,ierr)) 465 PetscCall(VecGetArrayF90(Top,top_v,ierr)) 466 PetscCall(VecGetArrayF90(Bottom,bottom_v,ierr)) 467 PetscCall(VecGetArrayF90(Left,left_v,ierr)) 468 PetscCall(VecGetArrayF90(Right,right_v,ierr)) 469 470! Initialize matrix entries to zero 471 PetscCall(MatAssembled(Hessian,assembled,ierr)) 472 if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(Hessian,ierr)) 473 474 rhx = real(mx) + 1.0 475 rhy = real(my) + 1.0 476 hx = 1.0/rhx 477 hy = 1.0/rhy 478 hydhx = hy/hx 479 hxdhy = hx/hy 480! compute Hessian over the locally owned part of the mesh 481 482 do i=xs,xs+xm-1 483 do j=ys,ys+ym-1 484 row = (j-gys)*gxm + (i-gxs) 485 486 xc = x_v(1+row) 487 xt = xc 488 xb = xc 489 xr = xc 490 xl = xc 491 xrb = xc 492 xlt = xc 493 494 if (i .eq. gxs) then ! Left side 495 xl = left_v(1+j - ys + 1) 496 xlt = left_v(1+j - ys + 2) 497 else 498 xl = x_v(1+row -1) 499 endif 500 501 if (j .eq. gys) then ! bottom side 502 xb = bottom_v(1+i - xs + 1) 503 xrb = bottom_v(1+i - xs + 2) 504 else 505 xb = x_v(1+row - gxm) 506 endif 507 508 if (i+1 .eq. gxs + gxm) then !right side 509 xr = right_v(1+j - ys + 1) 510 xrb = right_v(1+j - ys) 511 else 512 xr = x_v(1+row + 1) 513 endif 514 515 if (j+1 .eq. gym+gys) then !top side 516 xt = top_v(1+i - xs + 1) 517 xlt = top_v(1+i - xs) 518 else 519 xt = x_v(1+row + gxm) 520 endif 521 522 if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then 523 xlt = x_v(1+row - 1 + gxm) 524 endif 525 526 if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then 527 xrb = x_v(1+row + 1 - gxm) 528 endif 529 530 d1 = (xc-xl)*rhx 531 d2 = (xc-xr)*rhx 532 d3 = (xc-xt)*rhy 533 d4 = (xc-xb)*rhy 534 d5 = (xrb-xr)*rhy 535 d6 = (xrb-xb)*rhx 536 d7 = (xlt-xl)*rhy 537 d8 = (xlt-xt)*rhx 538 539 f1 = sqrt(1.0 + d1*d1 + d7*d7) 540 f2 = sqrt(1.0 + d1*d1 + d4*d4) 541 f3 = sqrt(1.0 + d3*d3 + d8*d8) 542 f4 = sqrt(1.0 + d3*d3 + d2*d2) 543 f5 = sqrt(1.0 + d2*d2 + d5*d5) 544 f6 = sqrt(1.0 + d4*d4 + d6*d6) 545 546 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2) 547 548 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4) 549 550 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4) 551 552 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2) 553 554 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6) 555 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3) 556 557 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) + & 558 & 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)+ & 559 & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4) 560 561 hl = hl * 0.5 562 hr = hr * 0.5 563 ht = ht * 0.5 564 hb = hb * 0.5 565 hbr = hbr * 0.5 566 htl = htl * 0.5 567 hc = hc * 0.5 568 569 k = 0 570 571 if (j .gt. 0) then 572 v(k) = hb 573 col(k) = row - gxm 574 k=k+1 575 endif 576 577 if ((j .gt. 0) .and. (i .lt. mx-1)) then 578 v(k) = hbr 579 col(k) = row-gxm+1 580 k=k+1 581 endif 582 583 if (i .gt. 0) then 584 v(k) = hl 585 col(k) = row - 1 586 k = k+1 587 endif 588 589 v(k) = hc 590 col(k) = row 591 k=k+1 592 593 if (i .lt. mx-1) then 594 v(k) = hr 595 col(k) = row + 1 596 k=k+1 597 endif 598 599 if ((i .gt. 0) .and. (j .lt. my-1)) then 600 v(k) = htl 601 col(k) = row + gxm - 1 602 k=k+1 603 endif 604 605 if (j .lt. my-1) then 606 v(k) = ht 607 col(k) = row + gxm 608 k=k+1 609 endif 610 611! Set matrix values using local numbering, defined earlier in main routine 612 PetscCall(MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES,ierr)) 613 614 enddo 615 enddo 616 617! restore vectors 618 PetscCall(VecRestoreArrayF90(localX,x_v,ierr)) 619 PetscCall(VecRestoreArrayF90(Left,left_v,ierr)) 620 PetscCall(VecRestoreArrayF90(Right,right_v,ierr)) 621 PetscCall(VecRestoreArrayF90(Top,top_v,ierr)) 622 PetscCall(VecRestoreArrayF90(Bottom,bottom_v,ierr)) 623 624! Assemble the matrix 625 PetscCall(MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr)) 626 PetscCall(MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr)) 627 628 PetscCall(PetscLogFlops(199.0d0*xm*ym,ierr)) 629 630 end 631 632! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h 633 634! ---------------------------------------------------------------------------- 635! 636!/* 637! MSA_BoundaryConditions - calculates the boundary conditions for the region 638! 639! 640!*/ 641 642 subroutine MSA_BoundaryConditions(ierr) 643 use plate2fmodule 644 implicit none 645 646 PetscErrorCode ierr 647 PetscInt i,j,k,limit,maxits 648 PetscInt xs, xm, gxs, gxm 649 PetscInt ys, ym, gys, gym 650 PetscInt bsize, lsize 651 PetscInt tsize, rsize 652 PetscReal one,two,three,tol 653 PetscReal scl,fnorm,det,xt 654 PetscReal yt,hx,hy,u1,u2,nf1,nf2 655 PetscReal njac11,njac12,njac21,njac22 656 PetscReal b, t, l, r 657 PetscReal,pointer :: boundary_v(:) 658 659 logical exitloop 660 PetscBool flg 661 662 limit=0 663 maxits = 5 664 tol=1e-10 665 b=-0.5 666 t= 0.5 667 l=-0.5 668 r= 0.5 669 xt=0 670 yt=0 671 one=1.0 672 two=2.0 673 three=3.0 674 675 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 676 PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)) 677 678 bsize = xm + 2 679 lsize = ym + 2 680 rsize = ym + 2 681 tsize = xm + 2 682 683 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr)) 684 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr)) 685 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr)) 686 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr)) 687 688 hx= (r-l)/(mx+1) 689 hy= (t-b)/(my+1) 690 691 do j=0,3 692 693 if (j.eq.0) then 694 yt=b 695 xt=l+hx*xs 696 limit=bsize 697 PetscCall(VecGetArrayF90(Bottom,boundary_v,ierr)) 698 699 elseif (j.eq.1) then 700 yt=t 701 xt=l+hx*xs 702 limit=tsize 703 PetscCall(VecGetArrayF90(Top,boundary_v,ierr)) 704 705 elseif (j.eq.2) then 706 yt=b+hy*ys 707 xt=l 708 limit=lsize 709 PetscCall(VecGetArrayF90(Left,boundary_v,ierr)) 710 711 elseif (j.eq.3) then 712 yt=b+hy*ys 713 xt=r 714 limit=rsize 715 PetscCall(VecGetArrayF90(Right,boundary_v,ierr)) 716 endif 717 718 do i=0,limit-1 719 720 u1=xt 721 u2=-yt 722 k = 0 723 exitloop = .false. 724 do while (k .lt. maxits .and. (.not. exitloop)) 725 726 nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt 727 nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt 728 fnorm=sqrt(nf1*nf1+nf2*nf2) 729 if (fnorm .gt. tol) then 730 njac11=one+u2*u2-u1*u1 731 njac12=two*u1*u2 732 njac21=-two*u1*u2 733 njac22=-one - u1*u1 + u2*u2 734 det = njac11*njac22-njac21*njac12 735 u1 = u1-(njac22*nf1-njac12*nf2)/det 736 u2 = u2-(njac11*nf2-njac21*nf1)/det 737 else 738 exitloop = .true. 739 endif 740 k=k+1 741 enddo 742 743 boundary_v(1+i) = u1*u1-u2*u2 744 if ((j .eq. 0) .or. (j .eq. 1)) then 745 xt = xt + hx 746 else 747 yt = yt + hy 748 endif 749 750 enddo 751 752 if (j.eq.0) then 753 PetscCall(VecRestoreArrayF90(Bottom,boundary_v,ierr)) 754 elseif (j.eq.1) then 755 PetscCall(VecRestoreArrayF90(Top,boundary_v,ierr)) 756 elseif (j.eq.2) then 757 PetscCall(VecRestoreArrayF90(Left,boundary_v,ierr)) 758 elseif (j.eq.3) then 759 PetscCall(VecRestoreArrayF90(Right,boundary_v,ierr)) 760 endif 761 762 enddo 763 764! Scale the boundary if desired 765 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-bottom',scl,flg,ierr)) 766 if (flg .eqv. PETSC_TRUE) then 767 PetscCall(VecScale(Bottom,scl,ierr)) 768 endif 769 770 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-top',scl,flg,ierr)) 771 if (flg .eqv. PETSC_TRUE) then 772 PetscCall(VecScale(Top,scl,ierr)) 773 endif 774 775 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-right',scl,flg,ierr)) 776 if (flg .eqv. PETSC_TRUE) then 777 PetscCall(VecScale(Right,scl,ierr)) 778 endif 779 780 PetscCall(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-left',scl,flg,ierr)) 781 if (flg .eqv. PETSC_TRUE) then 782 PetscCall(VecScale(Left,scl,ierr)) 783 endif 784 785 end 786 787! ---------------------------------------------------------------------------- 788! 789!/* 790! MSA_Plate - Calculates an obstacle for surface to stretch over 791! 792! Output Parameter: 793!. xl - lower bound vector 794!. xu - upper bound vector 795! 796!*/ 797 798 subroutine MSA_Plate(tao,xl,xu,dummy,ierr) 799 use plate2fmodule 800 implicit none 801 802 Tao tao 803 Vec xl,xu 804 PetscErrorCode ierr 805 PetscInt i,j,row 806 PetscInt xs, xm, ys, ym 807 PetscReal lb,ub 808 PetscInt dummy 809 PetscReal, pointer :: xl_v(:) 810 811 lb = PETSC_NINFINITY 812 ub = PETSC_INFINITY 813 814 if (bmy .lt. 0) bmy = 0 815 if (bmy .gt. my) bmy = my 816 if (bmx .lt. 0) bmx = 0 817 if (bmx .gt. mx) bmx = mx 818 819 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 820 821 PetscCall(VecSet(xl,lb,ierr)) 822 PetscCall(VecSet(xu,ub,ierr)) 823 824 PetscCall(VecGetArrayF90(xl,xl_v,ierr)) 825 826 do i=xs,xs+xm-1 827 828 do j=ys,ys+ym-1 829 830 row=(j-ys)*xm + (i-xs) 831 832 if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. & 833 & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then 834 xl_v(1+row) = bheight 835 836 endif 837 838 enddo 839 enddo 840 841 PetscCall(VecRestoreArrayF90(xl,xl_v,ierr)) 842 843 end 844 845! ---------------------------------------------------------------------------- 846! 847!/* 848! MSA_InitialPoint - Calculates an obstacle for surface to stretch over 849! 850! Output Parameter: 851!. X - vector for initial guess 852! 853!*/ 854 855 subroutine MSA_InitialPoint(X, ierr) 856 use plate2fmodule 857 implicit none 858 859 Vec X 860 PetscErrorCode ierr 861 PetscInt start,i,j 862 PetscInt row 863 PetscInt xs,xm,gxs,gxm 864 PetscInt ys,ym,gys,gym 865 PetscReal zero, np5 866 867 PetscReal,pointer :: left_v(:),right_v(:) 868 PetscReal,pointer :: bottom_v(:),top_v(:) 869 PetscReal,pointer :: x_v(:) 870 PetscBool flg 871 PetscRandom rctx 872 873 zero = 0.0 874 np5 = -0.5 875 876 PetscCall(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-start', start,flg,ierr)) 877 878 if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable 879 PetscCall(VecSet(X,zero,ierr)) 880 881 elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5 882 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)) 883 do i=0,start-1 884 PetscCall(VecSetRandom(X,rctx,ierr)) 885 enddo 886 887 PetscCall(PetscRandomDestroy(rctx,ierr)) 888 PetscCall(VecShift(X,np5,ierr)) 889 890 else ! average of boundary conditions 891 892! Get Local mesh boundaries 893 PetscCall(DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym,PETSC_NULL_INTEGER,ierr)) 894 PetscCall(DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,PETSC_NULL_INTEGER,ierr)) 895 896! Get pointers to vector data 897 PetscCall(VecGetArrayF90(Top,top_v,ierr)) 898 PetscCall(VecGetArrayF90(Bottom,bottom_v,ierr)) 899 PetscCall(VecGetArrayF90(Left,left_v,ierr)) 900 PetscCall(VecGetArrayF90(Right,right_v,ierr)) 901 902 PetscCall(VecGetArrayF90(localX,x_v,ierr)) 903 904! Perform local computations 905 do j=ys,ys+ym-1 906 do i=xs,xs+xm-1 907 row = (j-gys)*gxm + (i-gxs) 908 x_v(1+row) = ((j+1)*bottom_v(1+i-xs+1)/my + (my-j+1)*top_v(1+i-xs+1)/(my+2) + & 909 & (i+1)*left_v(1+j-ys+1)/mx + (mx-i+1)*right_v(1+j-ys+1)/(mx+2))*0.5 910 enddo 911 enddo 912 913! Restore vectors 914 PetscCall(VecRestoreArrayF90(localX,x_v,ierr)) 915 916 PetscCall(VecRestoreArrayF90(Left,left_v,ierr)) 917 PetscCall(VecRestoreArrayF90(Top,top_v,ierr)) 918 PetscCall(VecRestoreArrayF90(Bottom,bottom_v,ierr)) 919 PetscCall(VecRestoreArrayF90(Right,right_v,ierr)) 920 921 PetscCall(DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr)) 922 PetscCall(DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr)) 923 924 endif 925 926 end 927 928! 929!/*TEST 930! 931! build: 932! requires: !complex 933! 934! test: 935! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4 936! filter: sort -b 937! filter_output: sort -b 938! requires: !single 939! 940! test: 941! suffix: 2 942! nsize: 2 943! args: -tao_monitor_short -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4 944! filter: sort -b 945! filter_output: sort -b 946! requires: !single 947! 948!TEST*/ 949