1 /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */ 2 3 /* 4 Include "petsctao.h" so we can use TAO solvers. 5 petscdmda.h for distributed array 6 */ 7 #include <petsctao.h> 8 #include <petscdmda.h> 9 10 static char help[] = 11 "This example demonstrates use of the TAO package to \n\ 12 solve an unconstrained minimization problem. This example is based on a \n\ 13 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\ 14 boundary values along the edges of the domain, the objective is to find the\n\ 15 surface with the minimal area that satisfies the boundary conditions.\n\ 16 The command line options are:\n\ 17 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ 18 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ 19 -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\ 20 for an average of the boundary conditions\n\n"; 21 22 /* 23 User-defined application context - contains data needed by the 24 application-provided call-back routines, FormFunctionGradient() 25 and FormHessian(). 26 */ 27 typedef struct { 28 PetscInt mx, my; /* discretization in x, y directions */ 29 PetscReal *bottom, *top, *left, *right; /* boundary values */ 30 DM dm; /* distributed array data structure */ 31 Mat H; /* Hessian */ 32 } AppCtx; 33 34 /* -------- User-defined Routines --------- */ 35 36 static PetscErrorCode MSA_BoundaryConditions(AppCtx*); 37 static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec); 38 PetscErrorCode QuadraticH(AppCtx*,Vec,Mat); 39 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*); 40 PetscErrorCode FormGradient(Tao,Vec,Vec,void*); 41 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*); 42 PetscErrorCode My_Monitor(Tao, void *); 43 44 int main(int argc, char **argv) 45 { 46 PetscInt Nx, Ny; /* number of processors in x- and y- directions */ 47 Vec x; /* solution, gradient vectors */ 48 PetscBool flg, viewmat; /* flags */ 49 PetscBool fddefault, fdcoloring; /* flags */ 50 Tao tao; /* TAO solver context */ 51 AppCtx user; /* user-defined work context */ 52 ISColoring iscoloring; 53 MatFDColoring matfdcoloring; 54 55 /* Initialize TAO */ 56 PetscCall(PetscInitialize(&argc, &argv,(char *)0,help)); 57 58 /* Specify dimension of the problem */ 59 user.mx = 10; user.my = 10; 60 61 /* Check for any command line arguments that override defaults */ 62 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg)); 63 PetscCall(PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg)); 64 65 PetscCall(PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n")); 66 PetscCall(PetscPrintf(MPI_COMM_WORLD,"mx: %D my: %D \n\n",user.mx,user.my)); 67 68 /* Let PETSc determine the vector distribution */ 69 Nx = PETSC_DECIDE; Ny = PETSC_DECIDE; 70 71 /* Create distributed array (DM) to manage parallel grid and vectors */ 72 PetscCall(DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx, user.my,Nx,Ny,1,1,NULL,NULL,&user.dm)); 73 PetscCall(DMSetFromOptions(user.dm)); 74 PetscCall(DMSetUp(user.dm)); 75 76 /* Create TAO solver and set desired solution method.*/ 77 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 78 PetscCall(TaoSetType(tao,TAOCG)); 79 80 /* 81 Extract global vector from DA for the vector of variables -- PETSC routine 82 Compute the initial solution -- application specific, see below 83 Set this vector for use by TAO -- TAO routine 84 */ 85 PetscCall(DMCreateGlobalVector(user.dm,&x)); 86 PetscCall(MSA_BoundaryConditions(&user)); 87 PetscCall(MSA_InitialPoint(&user,x)); 88 PetscCall(TaoSetSolution(tao,x)); 89 90 /* 91 Initialize the Application context for use in function evaluations -- application specific, see below. 92 Set routines for function and gradient evaluation 93 */ 94 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 95 96 /* 97 Given the command line arguments, calculate the hessian with either the user- 98 provided function FormHessian, or the default finite-difference driven Hessian 99 functions 100 */ 101 PetscCall(PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault)); 102 PetscCall(PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring)); 103 104 /* 105 Create a matrix data structure to store the Hessian and set 106 the Hessian evalution routine. 107 Set the matrix structure to be used for Hessian evalutions 108 */ 109 PetscCall(DMCreateMatrix(user.dm,&user.H)); 110 PetscCall(MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE)); 111 112 if (fdcoloring) { 113 PetscCall(DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring)); 114 PetscCall(MatFDColoringCreate(user.H,iscoloring,&matfdcoloring)); 115 PetscCall(ISColoringDestroy(&iscoloring)); 116 PetscCall(MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user)); 117 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 118 PetscCall(TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring)); 119 } else if (fddefault) { 120 PetscCall(TaoSetHessian(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL)); 121 } else { 122 PetscCall(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user)); 123 } 124 125 /* 126 If my_monitor option is in command line, then use the user-provided 127 monitoring function 128 */ 129 PetscCall(PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat)); 130 if (viewmat) { 131 PetscCall(TaoSetMonitor(tao,My_Monitor,NULL,NULL)); 132 } 133 134 /* Check for any tao command line options */ 135 PetscCall(TaoSetFromOptions(tao)); 136 137 /* SOLVE THE APPLICATION */ 138 PetscCall(TaoSolve(tao)); 139 140 PetscCall(TaoView(tao,PETSC_VIEWER_STDOUT_WORLD)); 141 142 /* Free TAO data structures */ 143 PetscCall(TaoDestroy(&tao)); 144 145 /* Free PETSc data structures */ 146 PetscCall(VecDestroy(&x)); 147 PetscCall(MatDestroy(&user.H)); 148 if (fdcoloring) { 149 PetscCall(MatFDColoringDestroy(&matfdcoloring)); 150 } 151 PetscCall(PetscFree(user.bottom)); 152 PetscCall(PetscFree(user.top)); 153 PetscCall(PetscFree(user.left)); 154 PetscCall(PetscFree(user.right)); 155 PetscCall(DMDestroy(&user.dm)); 156 PetscCall(PetscFinalize()); 157 return 0; 158 } 159 160 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx) 161 { 162 PetscReal fcn; 163 164 PetscFunctionBegin; 165 PetscCall(FormFunctionGradient(tao,X,&fcn,G,userCtx)); 166 PetscFunctionReturn(0); 167 } 168 169 /* -------------------------------------------------------------------- */ 170 /* FormFunctionGradient - Evaluates the function and corresponding gradient. 171 172 Input Parameters: 173 . tao - the Tao context 174 . XX - input vector 175 . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient() 176 177 Output Parameters: 178 . fcn - the newly evaluated function 179 . GG - vector containing the newly evaluated gradient 180 */ 181 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx) 182 { 183 AppCtx *user = (AppCtx *) userCtx; 184 PetscInt i,j; 185 PetscInt mx=user->mx, my=user->my; 186 PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym; 187 PetscReal ft=0.0; 188 PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy; 189 PetscReal rhx=mx+1, rhy=my+1; 190 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 191 PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc; 192 PetscReal **g, **x; 193 Vec localX; 194 195 PetscFunctionBegin; 196 /* Get local mesh boundaries */ 197 PetscCall(DMGetLocalVector(user->dm,&localX)); 198 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 199 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 200 201 /* Scatter ghost points to local vector */ 202 PetscCall(DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX)); 203 PetscCall(DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX)); 204 205 /* Get pointers to vector data */ 206 PetscCall(DMDAVecGetArray(user->dm,localX,(void**)&x)); 207 PetscCall(DMDAVecGetArray(user->dm,G,(void**)&g)); 208 209 /* Compute function and gradient over the locally owned part of the mesh */ 210 for (j=ys; j<ys+ym; j++) { 211 for (i=xs; i< xs+xm; i++) { 212 213 xc = x[j][i]; 214 xlt=xrb=xl=xr=xb=xt=xc; 215 216 if (i==0) { /* left side */ 217 xl= user->left[j-ys+1]; 218 xlt = user->left[j-ys+2]; 219 } else { 220 xl = x[j][i-1]; 221 } 222 223 if (j==0) { /* bottom side */ 224 xb=user->bottom[i-xs+1]; 225 xrb =user->bottom[i-xs+2]; 226 } else { 227 xb = x[j-1][i]; 228 } 229 230 if (i+1 == gxs+gxm) { /* right side */ 231 xr=user->right[j-ys+1]; 232 xrb = user->right[j-ys]; 233 } else { 234 xr = x[j][i+1]; 235 } 236 237 if (j+1==gys+gym) { /* top side */ 238 xt=user->top[i-xs+1]; 239 xlt = user->top[i-xs]; 240 }else { 241 xt = x[j+1][i]; 242 } 243 244 if (i>gxs && j+1<gys+gym) { 245 xlt = x[j+1][i-1]; 246 } 247 if (j>gys && i+1<gxs+gxm) { 248 xrb = x[j-1][i+1]; 249 } 250 251 d1 = (xc-xl); 252 d2 = (xc-xr); 253 d3 = (xc-xt); 254 d4 = (xc-xb); 255 d5 = (xr-xrb); 256 d6 = (xrb-xb); 257 d7 = (xlt-xl); 258 d8 = (xt-xlt); 259 260 df1dxc = d1*hydhx; 261 df2dxc = (d1*hydhx + d4*hxdhy); 262 df3dxc = d3*hxdhy; 263 df4dxc = (d2*hydhx + d3*hxdhy); 264 df5dxc = d2*hydhx; 265 df6dxc = d4*hxdhy; 266 267 d1 *= rhx; 268 d2 *= rhx; 269 d3 *= rhy; 270 d4 *= rhy; 271 d5 *= rhy; 272 d6 *= rhx; 273 d7 *= rhy; 274 d8 *= rhx; 275 276 f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7); 277 f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4); 278 f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8); 279 f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2); 280 f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5); 281 f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6); 282 283 ft = ft + (f2 + f4); 284 285 df1dxc /= f1; 286 df2dxc /= f2; 287 df3dxc /= f3; 288 df4dxc /= f4; 289 df5dxc /= f5; 290 df6dxc /= f6; 291 292 g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5; 293 294 } 295 } 296 297 /* Compute triangular areas along the border of the domain. */ 298 if (xs==0) { /* left side */ 299 for (j=ys; j<ys+ym; j++) { 300 d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy; 301 d2=(user->left[j-ys+1] - x[j][0]) *rhx; 302 ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2); 303 } 304 } 305 if (ys==0) { /* bottom side */ 306 for (i=xs; i<xs+xm; i++) { 307 d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx; 308 d3=(user->bottom[i-xs+1]-x[0][i])*rhy; 309 ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2); 310 } 311 } 312 313 if (xs+xm==mx) { /* right side */ 314 for (j=ys; j< ys+ym; j++) { 315 d1=(x[j][mx-1] - user->right[j-ys+1])*rhx; 316 d4=(user->right[j-ys]-user->right[j-ys+1])*rhy; 317 ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4); 318 } 319 } 320 if (ys+ym==my) { /* top side */ 321 for (i=xs; i<xs+xm; i++) { 322 d1=(x[my-1][i] - user->top[i-xs+1])*rhy; 323 d4=(user->top[i-xs+1] - user->top[i-xs])*rhx; 324 ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4); 325 } 326 } 327 328 if (ys==0 && xs==0) { 329 d1=(user->left[0]-user->left[1])/hy; 330 d2=(user->bottom[0]-user->bottom[1])*rhx; 331 ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2); 332 } 333 if (ys+ym == my && xs+xm == mx) { 334 d1=(user->right[ym+1] - user->right[ym])*rhy; 335 d2=(user->top[xm+1] - user->top[xm])*rhx; 336 ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2); 337 } 338 339 ft=ft*area; 340 PetscCallMPI(MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD)); 341 342 /* Restore vectors */ 343 PetscCall(DMDAVecRestoreArray(user->dm,localX,(void**)&x)); 344 PetscCall(DMDAVecRestoreArray(user->dm,G,(void**)&g)); 345 346 /* Scatter values to global vector */ 347 PetscCall(DMRestoreLocalVector(user->dm,&localX)); 348 PetscCall(PetscLogFlops(67.0*xm*ym)); 349 PetscFunctionReturn(0); 350 } 351 352 /* ------------------------------------------------------------------- */ 353 /* 354 FormHessian - Evaluates Hessian matrix. 355 356 Input Parameters: 357 . tao - the Tao context 358 . x - input vector 359 . ptr - optional user-defined context, as set by TaoSetHessian() 360 361 Output Parameters: 362 . H - Hessian matrix 363 . Hpre - optionally different preconditioning matrix 364 . flg - flag indicating matrix structure 365 366 */ 367 PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr) 368 { 369 AppCtx *user = (AppCtx *) ptr; 370 371 PetscFunctionBegin; 372 /* Evaluate the Hessian entries*/ 373 PetscCall(QuadraticH(user,X,H)); 374 PetscFunctionReturn(0); 375 } 376 377 /* ------------------------------------------------------------------- */ 378 /* 379 QuadraticH - Evaluates Hessian matrix. 380 381 Input Parameters: 382 . user - user-defined context, as set by TaoSetHessian() 383 . X - input vector 384 385 Output Parameter: 386 . H - Hessian matrix 387 */ 388 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) 389 { 390 PetscInt i,j,k; 391 PetscInt mx=user->mx, my=user->my; 392 PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym; 393 PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy; 394 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 395 PetscReal hl,hr,ht,hb,hc,htl,hbr; 396 PetscReal **x, v[7]; 397 MatStencil col[7],row; 398 Vec localX; 399 PetscBool assembled; 400 401 PetscFunctionBegin; 402 /* Get local mesh boundaries */ 403 PetscCall(DMGetLocalVector(user->dm,&localX)); 404 405 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 406 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 407 408 /* Scatter ghost points to local vector */ 409 PetscCall(DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX)); 410 PetscCall(DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX)); 411 412 /* Get pointers to vector data */ 413 PetscCall(DMDAVecGetArray(user->dm,localX,(void**)&x)); 414 415 /* Initialize matrix entries to zero */ 416 PetscCall(MatAssembled(Hessian,&assembled)); 417 if (assembled) PetscCall(MatZeroEntries(Hessian)); 418 419 /* Set various matrix options */ 420 PetscCall(MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE)); 421 422 /* Compute Hessian over the locally owned part of the mesh */ 423 424 for (j=ys; j<ys+ym; j++) { 425 426 for (i=xs; i< xs+xm; i++) { 427 428 xc = x[j][i]; 429 xlt=xrb=xl=xr=xb=xt=xc; 430 431 /* Left side */ 432 if (i==0) { 433 xl = user->left[j-ys+1]; 434 xlt = user->left[j-ys+2]; 435 } else { 436 xl = x[j][i-1]; 437 } 438 439 if (j==0) { 440 xb = user->bottom[i-xs+1]; 441 xrb = user->bottom[i-xs+2]; 442 } else { 443 xb = x[j-1][i]; 444 } 445 446 if (i+1 == mx) { 447 xr = user->right[j-ys+1]; 448 xrb = user->right[j-ys]; 449 } else { 450 xr = x[j][i+1]; 451 } 452 453 if (j+1==my) { 454 xt = user->top[i-xs+1]; 455 xlt = user->top[i-xs]; 456 }else { 457 xt = x[j+1][i]; 458 } 459 460 if (i>0 && j+1<my) { 461 xlt = x[j+1][i-1]; 462 } 463 if (j>0 && i+1<mx) { 464 xrb = x[j-1][i+1]; 465 } 466 467 d1 = (xc-xl)/hx; 468 d2 = (xc-xr)/hx; 469 d3 = (xc-xt)/hy; 470 d4 = (xc-xb)/hy; 471 d5 = (xrb-xr)/hy; 472 d6 = (xrb-xb)/hx; 473 d7 = (xlt-xl)/hy; 474 d8 = (xlt-xt)/hx; 475 476 f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7); 477 f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4); 478 f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8); 479 f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2); 480 f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5); 481 f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6); 482 483 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ 484 (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2); 485 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ 486 (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4); 487 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ 488 (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4); 489 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ 490 (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2); 491 492 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6); 493 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3); 494 495 hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + 496 hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) + 497 (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) + 498 (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4); 499 500 hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0; 501 502 row.j = j; row.i = i; 503 k=0; 504 if (j>0) { 505 v[k]=hb; 506 col[k].j = j - 1; col[k].i = i; 507 k++; 508 } 509 510 if (j>0 && i < mx -1) { 511 v[k]=hbr; 512 col[k].j = j - 1; col[k].i = i+1; 513 k++; 514 } 515 516 if (i>0) { 517 v[k]= hl; 518 col[k].j = j; col[k].i = i-1; 519 k++; 520 } 521 522 v[k]= hc; 523 col[k].j = j; col[k].i = i; 524 k++; 525 526 if (i < mx-1) { 527 v[k]= hr; 528 col[k].j = j; col[k].i = i+1; 529 k++; 530 } 531 532 if (i>0 && j < my-1) { 533 v[k]= htl; 534 col[k].j = j+1; col[k].i = i-1; 535 k++; 536 } 537 538 if (j < my-1) { 539 v[k]= ht; 540 col[k].j = j+1; col[k].i = i; 541 k++; 542 } 543 544 /* 545 Set matrix values using local numbering, which was defined 546 earlier, in the main routine. 547 */ 548 PetscCall(MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES)); 549 } 550 } 551 552 PetscCall(DMDAVecRestoreArray(user->dm,localX,(void**)&x)); 553 PetscCall(DMRestoreLocalVector(user->dm,&localX)); 554 555 PetscCall(MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY)); 556 PetscCall(MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY)); 557 558 PetscCall(PetscLogFlops(199.0*xm*ym)); 559 PetscFunctionReturn(0); 560 } 561 562 /* ------------------------------------------------------------------- */ 563 /* 564 MSA_BoundaryConditions - Calculates the boundary conditions for 565 the region. 566 567 Input Parameter: 568 . user - user-defined application context 569 570 Output Parameter: 571 . user - user-defined application context 572 */ 573 static PetscErrorCode MSA_BoundaryConditions(AppCtx * user) 574 { 575 PetscInt i,j,k,limit=0,maxits=5; 576 PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym; 577 PetscInt mx=user->mx,my=user->my; 578 PetscInt bsize=0, lsize=0, tsize=0, rsize=0; 579 PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10; 580 PetscReal fnorm,det,hx,hy,xt=0,yt=0; 581 PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22; 582 PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5; 583 PetscReal *boundary; 584 PetscBool flg; 585 586 PetscFunctionBegin; 587 /* Get local mesh boundaries */ 588 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 589 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 590 591 bsize=xm+2; 592 lsize=ym+2; 593 rsize=ym+2; 594 tsize=xm+2; 595 596 PetscCall(PetscMalloc1(bsize,&user->bottom)); 597 PetscCall(PetscMalloc1(tsize,&user->top)); 598 PetscCall(PetscMalloc1(lsize,&user->left)); 599 PetscCall(PetscMalloc1(rsize,&user->right)); 600 601 hx= (r-l)/(mx+1); hy=(t-b)/(my+1); 602 603 for (j=0; j<4; j++) { 604 if (j==0) { 605 yt=b; 606 xt=l+hx*xs; 607 limit=bsize; 608 boundary=user->bottom; 609 } else if (j==1) { 610 yt=t; 611 xt=l+hx*xs; 612 limit=tsize; 613 boundary=user->top; 614 } else if (j==2) { 615 yt=b+hy*ys; 616 xt=l; 617 limit=lsize; 618 boundary=user->left; 619 } else { /* if (j==3) */ 620 yt=b+hy*ys; 621 xt=r; 622 limit=rsize; 623 boundary=user->right; 624 } 625 626 for (i=0; i<limit; i++) { 627 u1=xt; 628 u2=-yt; 629 for (k=0; k<maxits; k++) { 630 nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt; 631 nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt; 632 fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2); 633 if (fnorm <= tol) break; 634 njac11=one+u2*u2-u1*u1; 635 njac12=two*u1*u2; 636 njac21=-two*u1*u2; 637 njac22=-one - u1*u1 + u2*u2; 638 det = njac11*njac22-njac21*njac12; 639 u1 = u1-(njac22*nf1-njac12*nf2)/det; 640 u2 = u2-(njac11*nf2-njac21*nf1)/det; 641 } 642 643 boundary[i]=u1*u1-u2*u2; 644 if (j==0 || j==1) { 645 xt=xt+hx; 646 } else { /* if (j==2 || j==3) */ 647 yt=yt+hy; 648 } 649 650 } 651 652 } 653 654 /* Scale the boundary if desired */ 655 if (1==1) { 656 PetscReal scl = 1.0; 657 658 PetscCall(PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg)); 659 if (flg) { 660 for (i=0;i<bsize;i++) user->bottom[i]*=scl; 661 } 662 663 PetscCall(PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg)); 664 if (flg) { 665 for (i=0;i<tsize;i++) user->top[i]*=scl; 666 } 667 668 PetscCall(PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg)); 669 if (flg) { 670 for (i=0;i<rsize;i++) user->right[i]*=scl; 671 } 672 673 PetscCall(PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg)); 674 if (flg) { 675 for (i=0;i<lsize;i++) user->left[i]*=scl; 676 } 677 } 678 PetscFunctionReturn(0); 679 } 680 681 /* ------------------------------------------------------------------- */ 682 /* 683 MSA_InitialPoint - Calculates the initial guess in one of three ways. 684 685 Input Parameters: 686 . user - user-defined application context 687 . X - vector for initial guess 688 689 Output Parameters: 690 . X - newly computed initial guess 691 */ 692 static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X) 693 { 694 PetscInt start2=-1,i,j; 695 PetscReal start1=0; 696 PetscBool flg1,flg2; 697 698 PetscFunctionBegin; 699 PetscCall(PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1)); 700 PetscCall(PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2)); 701 702 if (flg1) { /* The zero vector is reasonable */ 703 704 PetscCall(VecSet(X, start1)); 705 706 } else if (flg2 && start2>0) { /* Try a random start between -0.5 and 0.5 */ 707 PetscRandom rctx; PetscReal np5=-0.5; 708 709 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rctx)); 710 PetscCall(VecSetRandom(X, rctx)); 711 PetscCall(PetscRandomDestroy(&rctx)); 712 PetscCall(VecShift(X, np5)); 713 714 } else { /* Take an average of the boundary conditions */ 715 PetscInt xs,xm,ys,ym; 716 PetscInt mx=user->mx,my=user->my; 717 PetscReal **x; 718 719 /* Get local mesh boundaries */ 720 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 721 722 /* Get pointers to vector data */ 723 PetscCall(DMDAVecGetArray(user->dm,X,(void**)&x)); 724 725 /* Perform local computations */ 726 for (j=ys; j<ys+ym; j++) { 727 for (i=xs; i< xs+xm; i++) { 728 x[j][i] = (((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0; 729 } 730 } 731 PetscCall(DMDAVecRestoreArray(user->dm,X,(void**)&x)); 732 PetscCall(PetscLogFlops(9.0*xm*ym)); 733 } 734 PetscFunctionReturn(0); 735 } 736 737 /*-----------------------------------------------------------------------*/ 738 PetscErrorCode My_Monitor(Tao tao, void *ctx) 739 { 740 Vec X; 741 742 PetscFunctionBegin; 743 PetscCall(TaoGetSolution(tao,&X)); 744 PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); 745 PetscFunctionReturn(0); 746 } 747 748 /*TEST 749 750 build: 751 requires: !complex 752 753 test: 754 args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3 755 requires: !single 756 757 test: 758 suffix: 2 759 nsize: 2 760 args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4 761 filter: grep -v "nls ksp" 762 requires: !single 763 764 test: 765 suffix: 3 766 nsize: 3 767 args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3 768 requires: !single 769 770 test: 771 suffix: 5 772 nsize: 2 773 args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3 774 requires: !single 775 776 TEST*/ 777