/*$Id: plate2.c 1.74 05/05/11 08:44:55-05:00 sarich@zorak.(none) $*/ #include "petscda.h" #include "tao.h" static char help[] = "This example demonstrates use of the TAO package to \n\ solve an unconstrained minimization problem. This example is based on a \n\ problem from the MINPACK-2 test suite. Given a rectangular 2-D domain, \n\ boundary values along the edges of the domain, and a plate represented by \n\ lower boundary conditions, the objective is to find the\n\ surface with the minimal area that satisfies the boundary conditions.\n\ The command line options are:\n\ -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\ -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\ -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\ -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\ -bheight <ht>, where <ht> = height of the plate\n\ -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\ for an average of the boundary conditions\n\n"; /*T Concepts: TAO - Solving a bound constrained minimization problem Routines: TaoInitialize(); TaoFinalize(); Routines: TaoApplicationCreate(); TaoAppDestroy(); Routines: TaoAppSetObjectiveAndGradientRoutine(); TaoAppSetHessianRoutine(); Routines: TaoAppSetInitialSolutionVec(); TaoAppSetHessianMat(); Routines: TaoAppSetVariableBounds(); Routines: TaoCreate(); TaoDestroy(); Routines: TaoSetOptions(); Processors: n T*/ /* User-defined application context - contains data needed by the application-provided call-back routines, FormFunctionGradient(), FormHessian(). */ typedef struct { /* problem parameters */ PetscReal bheight; /* Height of plate under the surface */ int mx, my; /* discretization in x, y directions */ int bmx,bmy; /* Size of plate under the surface */ Vec Bottom, Top, Left, Right; /* boundary values */ /* Working space */ Vec localX, localV; /* ghosted local vector */ DA da; /* distributed array data structure */ Mat H; } AppCtx; /* -------- User-defined Routines --------- */ static int MSA_BoundaryConditions(AppCtx*); static int MSA_InitialPoint(AppCtx*,Vec); static int MSA_Plate(TAO_APPLICATION,Vec,Vec,void*); int FormFunctionGradient(TAO_APPLICATION,Vec,double*,Vec,void*); int FormHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure*,void*); #undef __FUNCT__ #define __FUNCT__ "main" int main( int argc, char **argv ) { int info; /* used to check for functions returning nonzeros */ int Nx, Ny; /* number of processors in x- and y- directions */ int m, N; /* number of local and global elements in vectors */ Vec x; /* solution vector */ PetscTruth flg; /* A return variable when checking for user options */ TAO_SOLVER tao; /* TAO_SOLVER solver context */ TAO_APPLICATION plateapp; /* The PETSc application */ TaoMethod method = "tao_blmvm"; /* default minimization method */ double ff,gnorm,cnorm; /* iteration information */ int iter; ISLocalToGlobalMapping isltog; /* local-to-global mapping object */ TaoTerminateReason reason; AppCtx user; /* user-defined work context */ /* Initialize PETSc, TAO */ PetscInitialize( &argc, &argv,(char *)0,help ); TaoInitialize( &argc, &argv,(char *)0,help ); /* Specify default dimension of the problem */ user.mx = 10; user.my = 10; user.bheight=0.1; /* Check for any command line arguments that override defaults */ info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg); CHKERRQ(info); info = PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg); CHKERRQ(info); info = PetscOptionsGetReal(PETSC_NULL,"-bheight",&user.bheight,&flg); CHKERRQ(info); user.bmx = user.mx/2; user.bmy = user.my/2; info = PetscOptionsGetInt(PETSC_NULL,"-bmx",&user.bmx,&flg); CHKERRQ(info); info = PetscOptionsGetInt(PETSC_NULL,"-bmy",&user.bmy,&flg); CHKERRQ(info); PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n"); PetscPrintf(PETSC_COMM_WORLD,"mx:%d, my:%d, bmx:%d, bmy:%d, height:%4.2f\n", user.mx,user.my,user.bmx,user.bmy,user.bheight); /* Calculate any derived values from parameters */ N = user.mx*user.my; /* Let Petsc determine the dimensions of the local vectors */ Nx = PETSC_DECIDE; Ny = PETSC_DECIDE; /* A two dimensional distributed array will help define this problem, which derives from an elliptic PDE on two dimensional domain. From the distributed array, Create the vectors. */ info = DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx, user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da); CHKERRQ(info); /* Extract global and local vectors from DA; The local vectors are used solely as work space for the evaluation of the function, gradient, and Hessian. Duplicate for remaining vectors that are the same types. */ info = DACreateGlobalVector(user.da,&x); CHKERRQ(info); /* Solution */ info = DACreateLocalVector(user.da,&user.localX); CHKERRQ(info); info = VecDuplicate(user.localX,&user.localV); CHKERRQ(info); /* Create a matrix data structure to store the Hessian. Here we (optionally) also associate the local numbering scheme with the matrix so that later we can use local indices for matrix assembly. We could alternatively use global indices for matrix assembly. */ info = VecGetLocalSize(x,&m); CHKERRQ(info); info = MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL, 3,PETSC_NULL,&(user.H)); CHKERRQ(info); info = MatSetOption(user.H,MAT_SYMMETRIC); CHKERRQ(info); info = DAGetISLocalToGlobalMapping(user.da,&isltog); CHKERRQ(info); info = MatSetLocalToGlobalMapping(user.H,isltog); CHKERRQ(info); /* The TAO code begins here */ /* Create TAO solver and set desired solution method The method must either be 'tao_tron' or 'tao_blmvm' If blmvm is used, then hessian function is not called. */ info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info); info = TaoApplicationCreate(MPI_COMM_WORLD,&plateapp); CHKERRQ(info); /* Set initial solution guess; */ info = MSA_BoundaryConditions(&user); CHKERRQ(info); info = MSA_InitialPoint(&user,x); CHKERRQ(info); info = TaoAppSetInitialSolutionVec(plateapp,x); CHKERRQ(info); /* Set routines for function, gradient and hessian evaluation */ info = TaoAppSetObjectiveAndGradientRoutine(plateapp,FormFunctionGradient,(void*) &user); CHKERRQ(info); info = TaoAppSetHessianMat(plateapp,user.H,user.H); CHKERRQ(info); info = TaoAppSetHessianRoutine(plateapp,FormHessian,(void*)&user); CHKERRQ(info); /* Set Variable bounds */ info = TaoAppSetVariableBoundsRoutine(plateapp,MSA_Plate,(void*)&user); CHKERRQ(info); /* Check for any tao command line options */ info = TaoSetOptions(plateapp,tao); CHKERRQ(info); /* SOLVE THE APPLICATION */ info = TaoSolveApplication(plateapp,tao); CHKERRQ(info); /* Get information on termination */ info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,&cnorm,0,&reason); CHKERRQ(info); if (reason <= 0){ PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n"); PetscPrintf(MPI_COMM_WORLD,"Iteration: %d, f: %4.2e, Residual: %4.2e, Infeas: %4.2e\n",iter,ff,gnorm,cnorm); } /* Free TAO data structures */ info = TaoDestroy(tao); CHKERRQ(info); info = TaoAppDestroy(plateapp); CHKERRQ(info); /* Free PETSc data structures */ info = VecDestroy(x); CHKERRQ(info); info = MatDestroy(user.H); CHKERRQ(info); info = VecDestroy(user.localX); CHKERRQ(info); info = VecDestroy(user.localV); CHKERRQ(info); info = VecDestroy(user.Bottom); CHKERRQ(info); info = VecDestroy(user.Top); CHKERRQ(info); info = VecDestroy(user.Left); CHKERRQ(info); info = VecDestroy(user.Right); CHKERRQ(info); info = DADestroy(user.da); CHKERRQ(info); /* Finalize TAO and PETSc */ PetscFinalize(); TaoFinalize(); return 0; } #undef __FUNCT__ #define __FUNCT__ "FormFunctionGradient" /* FormFunctionGradient - Evaluates f(x) and gradient g(x). Input Parameters: . taoapp - the TAO_APPLICATION context . X - input vector . userCtx - optional user-defined context, as set by TaoSetPetscFunctionGradient() Output Parameters: . fcn - the function value . G - vector containing the newly evaluated gradient Notes: In this case, we discretize the domain and Create triangles. The surface of each triangle is planar, whose surface area can be easily computed. The total surface area is found by sweeping through the grid and computing the surface area of the two triangles that have their right angle at the grid point. The diagonal line segments on the grid that define the triangles run from top left to lower right. The numbering of points starts at the lower left and runs left to right, then bottom to top. */ int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn, Vec G,void *userCtx) { AppCtx * user = (AppCtx *) userCtx; int info,i,j,row; int mx=user->mx, my=user->my; int xs,xm,gxs,gxm,ys,ym,gys,gym; double ft=0; PetscReal zero=0.0; PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy; PetscReal rhx=mx+1, rhy=my+1; PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc; PetscReal *g, *x,*left,*right,*bottom,*top; Vec localX = user->localX, localG = user->localV; /* Get local mesh boundaries */ info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info); info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info); /* Scatter ghost points to local vector */ info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info); info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info); /* Initialize vector to zero */ info = VecSet(localG, zero); CHKERRQ(info); /* Get pointers to vector data */ info = VecGetArray(localX,&x); CHKERRQ(info); info = VecGetArray(localG,&g); CHKERRQ(info); info = VecGetArray(user->Top,&top); CHKERRQ(info); info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info); info = VecGetArray(user->Left,&left); CHKERRQ(info); info = VecGetArray(user->Right,&right); CHKERRQ(info); /* Compute function over the locally owned part of the mesh */ for (j=ys; j<ys+ym; j++){ for (i=xs; i< xs+xm; i++){ row=(j-gys)*gxm + (i-gxs); xc = x[row]; xlt=xrb=xl=xr=xb=xt=xc; if (i==0){ /* left side */ xl= left[j-ys+1]; xlt = left[j-ys+2]; } else { xl = x[row-1]; } if (j==0){ /* bottom side */ xb=bottom[i-xs+1]; xrb = bottom[i-xs+2]; } else { xb = x[row-gxm]; } if (i+1 == gxs+gxm){ /* right side */ xr=right[j-ys+1]; xrb = right[j-ys]; } else { xr = x[row+1]; } if (j+1==gys+gym){ /* top side */ xt=top[i-xs+1]; xlt = top[i-xs]; }else { xt = x[row+gxm]; } if (i>gxs && j+1<gys+gym){ xlt = x[row-1+gxm]; } if (j>gys && i+1<gxs+gxm){ xrb = x[row+1-gxm]; } d1 = (xc-xl); d2 = (xc-xr); d3 = (xc-xt); d4 = (xc-xb); d5 = (xr-xrb); d6 = (xrb-xb); d7 = (xlt-xl); d8 = (xt-xlt); df1dxc = d1*hydhx; df2dxc = ( d1*hydhx + d4*hxdhy ); df3dxc = d3*hxdhy; df4dxc = ( d2*hydhx + d3*hxdhy ); df5dxc = d2*hydhx; df6dxc = d4*hxdhy; d1 *= rhx; d2 *= rhx; d3 *= rhy; d4 *= rhy; d5 *= rhy; d6 *= rhx; d7 *= rhy; d8 *= rhx; f1 = sqrt( 1.0 + d1*d1 + d7*d7); f2 = sqrt( 1.0 + d1*d1 + d4*d4); f3 = sqrt( 1.0 + d3*d3 + d8*d8); f4 = sqrt( 1.0 + d3*d3 + d2*d2); f5 = sqrt( 1.0 + d2*d2 + d5*d5); f6 = sqrt( 1.0 + d4*d4 + d6*d6); ft = ft + (f2 + f4); df1dxc /= f1; df2dxc /= f2; df3dxc /= f3; df4dxc /= f4; df5dxc /= f5; df6dxc /= f6; g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5; } } /* Compute triangular areas along the border of the domain. */ if (xs==0){ /* left side */ for (j=ys; j<ys+ym; j++){ d3=(left[j-ys+1] - left[j-ys+2])*rhy; d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx; ft = ft+sqrt( 1.0 + d3*d3 + d2*d2); } } if (ys==0){ /* bottom side */ for (i=xs; i<xs+xm; i++){ d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx; d3=(bottom[i-xs+1]-x[i-gxs])*rhy; ft = ft+sqrt( 1.0 + d3*d3 + d2*d2); } } if (xs+xm==mx){ /* right side */ for (j=ys; j< ys+ym; j++){ d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx; d4=(right[j-ys]-right[j-ys+1])*rhy; ft = ft+sqrt( 1.0 + d1*d1 + d4*d4); } } if (ys+ym==my){ /* top side */ for (i=xs; i<xs+xm; i++){ d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy; d4=(top[i-xs+1] - top[i-xs])*rhx; ft = ft+sqrt( 1.0 + d1*d1 + d4*d4); } } if (ys==0 && xs==0){ d1=(left[0]-left[1])*rhy; d2=(bottom[0]-bottom[1])*rhx; ft +=sqrt( 1.0 + d1*d1 + d2*d2); } if (ys+ym == my && xs+xm == mx){ d1=(right[ym+1] - right[ym])*rhy; d2=(top[xm+1] - top[xm])*rhx; ft +=sqrt( 1.0 + d1*d1 + d2*d2); } ft=ft*area; info = MPI_Allreduce(&ft,fcn,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);CHKERRQ(info); /* Restore vectors */ info = VecRestoreArray(localX,&x); CHKERRQ(info); info = VecRestoreArray(localG,&g); CHKERRQ(info); info = VecRestoreArray(user->Left,&left); CHKERRQ(info); info = VecRestoreArray(user->Top,&top); CHKERRQ(info); info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info); info = VecRestoreArray(user->Right,&right); CHKERRQ(info); /* Scatter values to global vector */ info = DALocalToGlobal(user->da,localG,INSERT_VALUES,G); CHKERRQ(info); info = PetscLogFlops(70*xm*ym); CHKERRQ(info); return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "FormHessian" /* FormHessian - Evaluates Hessian matrix. Input Parameters: . taoapp - the TAO_APPLICATION context . x - input vector . ptr - optional user-defined context, as set by TaoSetHessian() Output Parameters: . A - Hessian matrix . B - optionally different preconditioning matrix . flag - flag indicating matrix structure Notes: Due to mesh point reordering with DAs, we must always work with the local mesh points, and then transform them to the new global numbering with the local-to-global mapping. We cannot work directly with the global numbers for the original uniprocessor mesh! Two methods are available for imposing this transformation when setting matrix entries: (A) MatSetValuesLocal(), using the local ordering (including ghost points!) - Do the following two steps once, before calling TaoSolve() - Use DAGetISLocalToGlobalMapping() to extract the local-to-global map from the DA - Associate this map with the matrix by calling MatSetLocalToGlobalMapping() - Then set matrix entries using the local ordering by calling MatSetValuesLocal() (B) MatSetValues(), using the global ordering - Use DAGetGlobalIndices() to extract the local-to-global map - Then apply this map explicitly yourself - Set matrix entries using the global ordering by calling MatSetValues() Option (A) seems cleaner/easier in many cases, and is the procedure used in this example. */ int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *Hptr, Mat *Hpc, MatStructure *flag, void *ptr) { int info; AppCtx *user = (AppCtx *) ptr; Mat Hessian = *Hpc; int i,j,k,row; int mx=user->mx, my=user->my; int xs,xm,gxs,gxm,ys,ym,gys,gym,col[7]; PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy; PetscReal rhx=mx+1, rhy=my+1; PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; PetscReal hl,hr,ht,hb,hc,htl,hbr; PetscReal *x,*left,*right,*bottom,*top; PetscReal v[7]; Vec localX = user->localX; PetscTruth assembled; /* Set various matrix options */ info = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES); CHKERRQ(info); info = MatSetOption(Hessian,MAT_COLUMNS_SORTED); CHKERRQ(info); info = MatSetOption(Hessian,MAT_ROWS_SORTED); CHKERRQ(info); /* Initialize matrix entries to zero */ info = MatAssembled(Hessian,&assembled); CHKERRQ(info); if (assembled){info = MatZeroEntries(Hessian); CHKERRQ(info);} /* Get local mesh boundaries */ info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info); info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info); /* Scatter ghost points to local vector */ info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info); info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info); /* Get pointers to vector data */ info = VecGetArray(localX,&x); CHKERRQ(info); info = VecGetArray(user->Top,&top); CHKERRQ(info); info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info); info = VecGetArray(user->Left,&left); CHKERRQ(info); info = VecGetArray(user->Right,&right); CHKERRQ(info); /* Compute Hessian over the locally owned part of the mesh */ for (i=xs; i< xs+xm; i++){ for (j=ys; j<ys+ym; j++){ row=(j-gys)*gxm + (i-gxs); xc = x[row]; xlt=xrb=xl=xr=xb=xt=xc; /* Left side */ if (i==gxs){ xl= left[j-ys+1]; xlt = left[j-ys+2]; } else { xl = x[row-1]; } if (j==gys){ xb=bottom[i-xs+1]; xrb = bottom[i-xs+2]; } else { xb = x[row-gxm]; } if (i+1 == gxs+gxm){ xr=right[j-ys+1]; xrb = right[j-ys]; } else { xr = x[row+1]; } if (j+1==gys+gym){ xt=top[i-xs+1]; xlt = top[i-xs]; }else { xt = x[row+gxm]; } if (i>gxs && j+1<gys+gym){ xlt = x[row-1+gxm]; } if (j>gys && i+1<gxs+gxm){ xrb = x[row+1-gxm]; } d1 = (xc-xl)*rhx; d2 = (xc-xr)*rhx; d3 = (xc-xt)*rhy; d4 = (xc-xb)*rhy; d5 = (xrb-xr)*rhy; d6 = (xrb-xb)*rhx; d7 = (xlt-xl)*rhy; d8 = (xlt-xt)*rhx; f1 = sqrt( 1.0 + d1*d1 + d7*d7); f2 = sqrt( 1.0 + d1*d1 + d4*d4); f3 = sqrt( 1.0 + d3*d3 + d8*d8); f4 = sqrt( 1.0 + d3*d3 + d2*d2); f5 = sqrt( 1.0 + d2*d2 + d5*d5); f6 = sqrt( 1.0 + d4*d4 + d6*d6); hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2); hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4); ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4); hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2); hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6); htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3); 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) + 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)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4); hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5; k=0; if (j>0){ v[k]=hb; col[k]=row - gxm; k++; } if (j>0 && i < mx -1){ v[k]=hbr; col[k]=row - gxm+1; k++; } if (i>0){ v[k]= hl; col[k]=row - 1; k++; } v[k]= hc; col[k]=row; k++; if (i < mx-1 ){ v[k]= hr; col[k]=row+1; k++; } if (i>0 && j < my-1 ){ v[k]= htl; col[k] = row+gxm-1; k++; } if (j < my-1 ){ v[k]= ht; col[k] = row+gxm; k++; } /* Set matrix values using local numbering, which was defined earlier, in the main routine. */ info = MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info); } } /* Restore vectors */ info = VecRestoreArray(localX,&x); CHKERRQ(info); info = VecRestoreArray(user->Left,&left); CHKERRQ(info); info = VecRestoreArray(user->Top,&top); CHKERRQ(info); info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info); info = VecRestoreArray(user->Right,&right); CHKERRQ(info); /* Assemble the matrix */ info = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info); info = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info); info = PetscLogFlops(199*xm*ym); CHKERRQ(info); return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "MSA_BoundaryConditions" /* MSA_BoundaryConditions - Calculates the boundary conditions for the region. Input Parameter: . user - user-defined application context Output Parameter: . user - user-defined application context */ static int MSA_BoundaryConditions(AppCtx * user) { int i,j,k,limit=0,info,maxits=5; int xs,ys,xm,ym,gxs,gys,gxm,gym; int mx=user->mx,my=user->my; int bsize=0, lsize=0, tsize=0, rsize=0; PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10; PetscReal fnorm,det,hx,hy,xt=0,yt=0; PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22; PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5; PetscReal *boundary; PetscTruth flg; Vec Bottom,Top,Right,Left; /* Get local mesh boundaries */ info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info); info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info); bsize=xm+2; lsize=ym+2; rsize=ym+2; tsize=xm+2; info = VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom); CHKERRQ(info); info = VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top); CHKERRQ(info); info = VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left); CHKERRQ(info); info = VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right); CHKERRQ(info); user->Top=Top; user->Left=Left; user->Bottom=Bottom; user->Right=Right; hx= (r-l)/(mx+1); hy=(t-b)/(my+1); for (j=0; j<4; j++){ if (j==0){ yt=b; xt=l+hx*xs; limit=bsize; VecGetArray(Bottom,&boundary); } else if (j==1){ yt=t; xt=l+hx*xs; limit=tsize; VecGetArray(Top,&boundary); } else if (j==2){ yt=b+hy*ys; xt=l; limit=lsize; VecGetArray(Left,&boundary); } else if (j==3){ yt=b+hy*ys; xt=r; limit=rsize; VecGetArray(Right,&boundary); } for (i=0; i<limit; i++){ u1=xt; u2=-yt; for (k=0; k<maxits; k++){ nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt; nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt; fnorm=sqrt(nf1*nf1+nf2*nf2); if (fnorm <= tol) break; njac11=one+u2*u2-u1*u1; njac12=two*u1*u2; njac21=-two*u1*u2; njac22=-one - u1*u1 + u2*u2; det = njac11*njac22-njac21*njac12; u1 = u1-(njac22*nf1-njac12*nf2)/det; u2 = u2-(njac11*nf2-njac21*nf1)/det; } boundary[i]=u1*u1-u2*u2; if (j==0 || j==1) { xt=xt+hx; } else if (j==2 || j==3){ yt=yt+hy; } } if (j==0){ info = VecRestoreArray(Bottom,&boundary); CHKERRQ(info); } else if (j==1){ info = VecRestoreArray(Top,&boundary); CHKERRQ(info); } else if (j==2){ info = VecRestoreArray(Left,&boundary); CHKERRQ(info); } else if (j==3){ info = VecRestoreArray(Right,&boundary); CHKERRQ(info); } } /* Scale the boundary if desired */ info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg); CHKERRQ(info); if (flg){ info = VecScale(Bottom, scl); CHKERRQ(info); } info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg); CHKERRQ(info); if (flg){ info = VecScale(Top, scl); CHKERRQ(info); } info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg); CHKERRQ(info); if (flg){ info = VecScale(Right, scl); CHKERRQ(info); } info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg); CHKERRQ(info); if (flg){ info = VecScale(Left, scl); CHKERRQ(info); } return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "MSA_Plate" /* MSA_Plate - Calculates an obstacle for surface to stretch over. Input Parameter: . user - user-defined application context Output Parameter: . user - user-defined application context */ static int MSA_Plate(TAO_APPLICATION taoapp, Vec XL,Vec XU,void *ctx){ AppCtx *user=(AppCtx *)ctx; int i,j,row,info; int xs,ys,xm,ym; int mx=user->mx, my=user->my, bmy, bmx; PetscReal t1,t2,t3; PetscReal *xl, lb=TAO_NINFINITY, ub=TAO_INFINITY; PetscTruth cylinder; user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy); user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx); bmy=user->bmy, bmx=user->bmx; info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info); info = VecSet(XL, lb); CHKERRQ(info); info = VecSet(XU, ub); CHKERRQ(info); info = VecGetArray(XL,&xl); CHKERRQ(info); info = PetscOptionsHasName(PETSC_NULL,"-cylinder",&cylinder); CHKERRQ(info); /* Compute the optional lower box */ if (cylinder){ for (i=xs; i< xs+xm; i++){ for (j=ys; j<ys+ym; j++){ row=(j-ys)*xm + (i-xs); t1=(2.0*i-mx)*bmy; t2=(2.0*j-my)*bmx; t3=bmx*bmx*bmy*bmy; if ( t1*t1 + t2*t2 <= t3 ){ xl[row] = user->bheight; } } } } else { /* Compute the optional lower box */ for (i=xs; i< xs+xm; i++){ for (j=ys; j<ys+ym; j++){ row=(j-ys)*xm + (i-xs); if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 && j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){ xl[row] = user->bheight; } } } } info = VecRestoreArray(XL,&xl); CHKERRQ(info); return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "MSA_InitialPoint" /* MSA_InitialPoint - Calculates the initial guess in one of three ways. Input Parameters: . user - user-defined application context . X - vector for initial guess Output Parameters: . X - newly computed initial guess */ static int MSA_InitialPoint(AppCtx * user, Vec X) { int start=-1,i,j,info; PetscReal zero=0.0; PetscTruth flg; info = PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); CHKERRQ(info); if (flg && start==0){ /* The zero vector is reasonable */ info = VecSet(X, zero); CHKERRQ(info); /* PetscLogInfo((user,"Min. Surface Area Problem: Start with 0 vector \n")); */ } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */ PetscRandom rctx; PetscReal np5=-0.5; info = PetscRandomCreate(MPI_COMM_WORLD,RANDOM_DEFAULT,&rctx); CHKERRQ(info); for (i=0; i<start; i++){ info = VecSetRandom(X, rctx); CHKERRQ(info); } info = PetscRandomDestroy(rctx); CHKERRQ(info); info = VecShift(X, np5); /* PetscLogInfo((user,"Min. Surface Area Problem: Start with random vector %d\n",start)); */ } else { /* Take an average of the boundary conditions */ int row,xs,xm,gxs,gxm,ys,ym,gys,gym; int mx=user->mx,my=user->my; PetscReal *x,*left,*right,*bottom,*top; Vec localX = user->localX; /* Get local mesh boundaries */ info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info); info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info); /* Get pointers to vector data */ info = VecGetArray(user->Top,&top); CHKERRQ(info); info = VecGetArray(user->Bottom,&bottom); CHKERRQ(info); info = VecGetArray(user->Left,&left); CHKERRQ(info); info = VecGetArray(user->Right,&right); CHKERRQ(info); info = VecGetArray(localX,&x); CHKERRQ(info); /* Perform local computations */ for (j=ys; j<ys+ym; j++){ for (i=xs; i< xs+xm; i++){ row=(j-gys)*gxm + (i-gxs); x[row] = ( (j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+ (i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0; } } /* Restore vectors */ info = VecRestoreArray(localX,&x); CHKERRQ(info); info = VecRestoreArray(user->Left,&left); CHKERRQ(info); info = VecRestoreArray(user->Top,&top); CHKERRQ(info); info = VecRestoreArray(user->Bottom,&bottom); CHKERRQ(info); info = VecRestoreArray(user->Right,&right); CHKERRQ(info); /* Scatter values into global vector */ info = DALocalToGlobal(user->da,localX,INSERT_VALUES,X); CHKERRQ(info); } return 0; }