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