! "$Id: plate2f.F 1.47 05/05/12 13:59:03-05:00 sarich@zorak.(none) $"; ! ! Program usage: mpirun -np <proc> plate2f [all TAO options] ! ! This example demonstrates use of the TAO package to solve a bound constrained ! minimization problem. This example is based on a problem from the ! MINPACK-2 test suite. Given a rectangular 2-D domain and boundary values ! along the edges of the domain, the objective is to find the surface ! with the minimal area that satisfies the boundary conditions. ! The command line options are: ! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction ! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction ! -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction ! -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction ! -bheight <ht>, where <ht> = height of the plate ! !/*T ! Concepts: TAO - Solving a bound constrained minimization problem ! Routines: TaoInitialize(); TaoFinalize(); ! Routines: TaoCreate(); TaoDestroy(); ! Routines: TaoAppSetObjectiveAndGradientRoutine(); ! Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine(); ! Routines: TaoAppSetInitialSolutionVec(); TaoAppSetVariableBounds(); ! Routines: TaoSetApplication(); TaoSetOptions(); ! Routines: TaoApplicationCreate(); TaoSolve(); ! Routines: TaoView(); TaoAppDestroy(); ! Processors: n !T*/ implicit none #include "plate2f.h" ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Variable declarations ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! Variables: ! (common from plate2f.h): ! Nx, Ny number of processors in x- and y- directions ! mx, my number of grid points in x,y directions ! N global dimension of vector integer info ! used to check for functions returning nonzeros Vec x ! solution vector Vec xl, xu ! lower and upper bounds vectorsp integer m ! number of local elements in vector TAO_SOLVER tao ! TAO_SOLVER solver context TAO_APPLICATION plateapp ! PETSc application Mat H ! Hessian matrix ISLocalToGlobalMapping isltog ! local to global mapping object PetscTruth flg external FormFunctionGradient external FormHessian external MSA_BoundaryConditions external MSA_Plate external MSA_InitialPoint ! Initialize Tao call PetscInitialize(PETSC_NULL_CHARACTER,info) call TaoInitialize(PETSC_NULL_CHARACTER,info) ! Specify default dimensions of the problem mx = 10 my = 10 bheight = 0.1 ! Check for any command line arguments that override defaults call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-mx",mx,flg,info) call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-my",my,flg,info) bmx = mx/2 bmy = my/2 call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmx",bmx,flg,info) call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmy",bmy,flg,info) call PetscOptionsGetReal(TAO_NULL_CHARACTER,"-bheight",bheight, & & flg,info) ! Calculate any derived values from parameters N = mx*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 a two-dimensional domain. From the ! distributed array, create the vectors call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX, & & mx,my,Nx,Ny,1,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & & da,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. call DACreateGlobalVector(da,x,info) call DACreateLocalVector(da,localX,info) call VecDuplicate(localX,localV,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 call VecGetLocalSize(x,m,info) call MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL_INTEGER, & & 3,PETSC_NULL_INTEGER,H,info) call MatSetOption(H,MAT_SYMMETRIC,info) call DAGetISLocalToGlobalMapping(da,isltog,info) call MatSetLocalToGlobalMapping(H,isltog,info) ! The Tao code begins here ! Create TAO solver and set desired solution method. ! This problems uses bounded variables, so the ! method must either be 'tao_tron' or 'tao_blmvm' call TaoCreate(MPI_COMM_WORLD,'tao_blmvm',tao,info) call TaoApplicationCreate(MPI_COMM_WORLD,plateapp,info) ! Set minimization function and gradient, hessian evaluation functions ! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers call TaoAppSetObjectiveAndGradientRo(plateapp, & & FormFunctionGradient,PETSC_NULL_OBJECT,info) call TaoAppSetHessianMat(plateapp,H,H,info) call TaoAppSetHessianRoutine(plateapp,FormHessian, & & PETSC_NULL_OBJECT, info) ! Set Variable bounds call MSA_BoundaryConditions(info) call VecDuplicate(x,xl,info) call VecDuplicate(x,xu,info) call MSA_Plate(xl,xu,info) call TaoAppSetVariableBounds(plateapp,xl,xu,info) ! Set the initial solution guess call MSA_InitialPoint(x, info) call TaoAppSetInitialSolutionVec(plateapp,x,info) ! Now that the PETSc application is all set, attach to tao context call TaoSetApplication(tao,plateapp,info) ! Check for any tao command line options call TaoSetOptions(plateapp,tao,info) ! Solve the application call TaoSolveApplication(plateapp,tao,info) ! View TAO solver information ! call TaoView(tao,info) ! Free TAO data structures call TaoDestroy(tao,info) call TaoAppDestroy(plateapp,info) ! Free PETSc data structures call VecDestroy(x,info) call VecDestroy(xl,info) call VecDestroy(xu,info) call VecDestroy(Top,info) call VecDestroy(Bottom,info) call VecDestroy(Left,info) call VecDestroy(Right,info) call MatDestroy(H,info) call VecDestroy(localX,info) call VecDestroy(localV,info) call DADestroy(da,info) ! Finalize TAO call TaoFinalize(info) call PetscFinalize(info) end ! --------------------------------------------------------------------- ! ! FormFunctionGradient - Evaluates function f(X). ! ! Input Parameters: ! tao - the TAO_SOLVER context ! X - the input vector ! dummy - optional user-defined context, as set by TaoSetFunction() ! (not used here) ! ! Output Parameters: ! fcn - the newly evaluated function ! G - the gradient vector ! info - error code ! subroutine FormFunctionGradient(tao,X,fcn,G,dummy,info) implicit none ! da, localX, localG, Top, Bottom, Left, Right defined in plate2f.h #include "plate2f.h" ! Input/output variables TAO_SOLVER tao PetscScalar fcn Vec X, G integer dummy, info integer i,j,row integer xs, xm, gxs, gxm, ys, ym, gys, gym PetscScalar ft,zero,hx,hy,hydhx,hxdhy,area,rhx,rhy PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8 PetscScalar df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc PetscScalar xc,xl,xr,xt,xb,xlt,xrb ! PETSc's VecGetArray acts differently in Fortran than it does in C. ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info) ! will return an array of doubles referenced by x_array offset by x_index. ! i.e., to reference the kth element of X, use x_array(k + x_index). ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. PetscScalar g_v(0:1),x_v(0:1), top_v(0:1) PetscScalar left_v(0:1), right_v(0:1),bottom_v(0:1) PetscOffset g_i,left_i,right_i,bottom_i,top_i PetscOffset x_i ft = 0.0d0 zero = 0.0d0 hx = 1.0d0/(mx + 1) hy = 1.0d0/(my + 1) hydhx = hy/hx hxdhy = hx/hy area = 0.5d0 * hx * hy rhx = mx + 1.0d0 rhy = my + 1.0d0 ! Get local mesh boundaries call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, & & PETSC_NULL_INTEGER,info) call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, & & gxm,gym,PETSC_NULL_INTEGER,info) ! Scatter ghost points to local vector call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info) call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info) ! Initialize the vector to zero call VecSet(localV,zero,info) ! Get arrays to vector data (See note above about using VecGetArray in Fortran) call VecGetArray(localX,x_v,x_i,info) call VecGetArray(localV,g_v,g_i,info) call VecGetArray(Top,top_v,top_i,info) call VecGetArray(Bottom,bottom_v,bottom_i,info) call VecGetArray(Left,left_v,left_i,info) call VecGetArray(Right,right_v,right_i,info) ! Compute function over the locally owned part of the mesh do j = ys,ys+ym-1 do i = xs,xs+xm-1 row = (j-gys)*gxm + (i-gxs) xc = x_v(row+x_i) xt = xc xb = xc xr = xc xl = xc xrb = xc xlt = xc if (i .eq. 0) then !left side xl = left_v(j - ys + 1 + left_i) xlt = left_v(j - ys + 2 + left_i) else xl = x_v(row - 1 + x_i) endif if (j .eq. 0) then !bottom side xb = bottom_v(i - xs + 1 + bottom_i) xrb = bottom_v(i - xs + 2 + bottom_i) else xb = x_v(row - gxm + x_i) endif if (i + 1 .eq. gxs + gxm) then !right side xr = right_v(j - ys + 1 + right_i) xrb = right_v(j - ys + right_i) else xr = x_v(row + 1 + x_i) endif if (j + 1 .eq. gys + gym) then !top side xt = top_v(i - xs + 1 + top_i) xlt = top_v(i - xs + top_i) else xt = x_v(row + gxm + x_i) endif if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then xlt = x_v(row - 1 + gxm + x_i) endif if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then xrb = x_v(row + 1 - gxm + x_i) endif 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 = d1 * rhx d2 = d2 * rhx d3 = d3 * rhy d4 = d4 * rhy d5 = d5 * rhy d6 = d6 * rhx d7 = d7 * rhy d8 = d8 * rhx f1 = sqrt(1.0d0 + d1*d1 + d7*d7) f2 = sqrt(1.0d0 + d1*d1 + d4*d4) f3 = sqrt(1.0d0 + d3*d3 + d8*d8) f4 = sqrt(1.0d0 + d3*d3 + d2*d2) f5 = sqrt(1.0d0 + d2*d2 + d5*d5) f6 = sqrt(1.0d0 + d4*d4 + d6*d6) ft = ft + f2 + f4 df1dxc = df1dxc / f1 df2dxc = df2dxc / f2 df3dxc = df3dxc / f3 df4dxc = df4dxc / f4 df5dxc = df5dxc / f5 df6dxc = df6dxc / f6 g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + & & df5dxc + df6dxc) enddo enddo ! Compute triangular areas along the border of the domain. if (xs .eq. 0) then ! left side do j=ys,ys+ym-1 d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) & & * rhy d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) & & * rhx ft = ft + sqrt(1.0d0 + d3*d3 + d2*d2) enddo endif if (ys .eq. 0) then !bottom side do i=xs,xs+xm-1 d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) & & * rhx d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy ft = ft + sqrt(1.0 + d3*d3 + d2*d2) enddo endif if (xs + xm .eq. mx) then ! right side do j=ys,ys+ym-1 d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4) enddo endif if (ys + ym .eq. my) then do i=xs,xs+xm-1 d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4) enddo endif if ((ys .eq. 0) .and. (xs .eq. 0)) then d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2) endif if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2) endif ft = ft * area call MPI_Allreduce(ft,fcn,1,MPI_DOUBLE_PRECISION, & & MPI_SUM,MPI_COMM_WORLD,info) ! Restore vectors call VecRestoreArray(localX,x_v,x_i,info) call VecRestoreArray(localV,g_v,g_i,info) call VecRestoreArray(Left,left_v,left_i,info) call VecRestoreArray(Top,top_v,top_i,info) call VecRestoreArray(Bottom,bottom_v,bottom_i,info) call VecRestoreArray(Right,right_v,right_i,info) ! Scatter values to global vector call DALocalToGlobal(da,localV,INSERT_VALUES,G,info) call PetscLogFlops(70*xm*ym,info) return end !FormFunctionGradient ! ---------------------------------------------------------------------------- ! !/* ! FormHessian - Evaluates Hessian matrix. ! ! Input Parameters: !. tao - the TAO_SOLVER context !. X - input vector !. dummy - not used ! ! Output Parameters: !. Hessian - Hessian matrix !. Hpc - 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. */ subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, info) implicit none ! da,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h #include "plate2f.h" TAO_SOLVER tao Vec X Mat Hessian,Hpc MatStructure flg integer dummy,info integer i,j,k,row integer xs,xm,gxs,gxm,ys,ym,gys,gym,col(0:6) PetscScalar hx,hy,hydhx,hxdhy,rhx,rhy PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8 PetscScalar xc,xl,xr,xt,xb,xlt,xrb PetscScalar hl,hr,ht,hb,hc,htl,hbr ! PETSc's VecGetArray acts differently in Fortran than it does in C. ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info) ! will return an array of doubles referenced by x_array offset by x_index. ! i.e., to reference the kth element of X, use x_array(k + x_index). ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. PetscScalar right_v(0:1),left_v(0:1),bottom_v(0:1),top_v(0:1) PetscScalar x_v(0:1) PetscOffset x_i,right_i,left_i,bottom_i,top_i PetscScalar v(0:6) PetscTruth assembled ! Set various matrix options call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,info) call MatSetOption(Hessian,MAT_COLUMNS_SORTED,info) call MatSetOption(Hessian,MAT_ROWS_SORTED,info) ! Get local mesh boundaries call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, & & PETSC_NULL_INTEGER,info) call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, & & PETSC_NULL_INTEGER,info) ! Scatter ghost points to local vectors call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info) call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info) ! Get pointers to vector data (see note on Fortran arrays above) call VecGetArray(localX,x_v,x_i,info) call VecGetArray(Top,top_v,top_i,info) call VecGetArray(Bottom,bottom_v,bottom_i,info) call VecGetArray(Left,left_v,left_i,info) call VecGetArray(Right,right_v,right_i,info) ! Initialize matrix entries to zero call MatAssembled(Hessian,assembled,info) if (assembled .eq. PETSC_TRUE) call MatZeroEntries(Hessian,info) rhx = mx + 1.0 rhy = my + 1.0 hx = 1.0/rhx hy = 1.0/rhy hydhx = hy/hx hxdhy = hx/hy ! compute Hessian over the locally owned part of the mesh do i=xs,xs+xm-1 do j=ys,ys+ym-1 row = (j-gys)*gxm + (i-gxs) xc = x_v(row + x_i) xt = xc xb = xc xr = xc xl = xc xrb = xc xlt = xc if (i .eq. gxs) then ! Left side xl = left_v(left_i + j - ys + 1) xlt = left_v(left_i + j - ys + 2) else xl = x_v(x_i + row -1 ) endif if (j .eq. gys) then ! bottom side xb = bottom_v(bottom_i + i - xs + 1) xrb = bottom_v(bottom_i + i - xs + 2) else xb = x_v(x_i + row - gxm) endif if (i+1 .eq. gxs + gxm) then !right side xr = right_v(right_i + j - ys + 1) xrb = right_v(right_i + j - ys) else xr = x_v(x_i + row + 1) endif if (j+1 .eq. gym+gys) then !top side xt = top_v(top_i +i - xs + 1) xlt = top_v(top_i + i - xs) else xt = x_v(x_i + row + gxm) endif if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then xlt = x_v(x_i + row - 1 + gxm) endif if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then xrb = x_v(x_i + row + 1 - gxm) endif 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.0d0 + d1*d1 + d7*d7) f2 = sqrt( 1.0d0 + d1*d1 + d4*d4) f3 = sqrt( 1.0d0 + d3*d3 + d8*d8) f4 = sqrt( 1.0d0 + d3*d3 + d2*d2) f5 = sqrt( 1.0d0 + d2*d2 + d5*d5) f6 = sqrt( 1.0d0 + 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 = hl * 0.5 hr = hr * 0.5 ht = ht * 0.5 hb = hb * 0.5 hbr = hbr * 0.5 htl = htl * 0.5 hc = hc * 0.5 k = 0 if (j .gt. 0) then v(k) = hb col(k) = row - gxm k=k+1 endif if ((j .gt. 0) .and. (i .lt. mx-1)) then v(k) = hbr col(k) = row-gxm+1 k=k+1 endif if (i .gt. 0) then v(k) = hl col(k) = row - 1 k = k+1 endif v(k) = hc col(k) = row k=k+1 if (i .lt. mx-1) then v(k) = hr col(k) = row + 1 k=k+1 endif if ((i .gt. 0) .and. (j .lt. my-1)) then v(k) = htl col(k) = row + gxm - 1 k=k+1 endif if (j .lt. my-1) then v(k) = ht col(k) = row + gxm k=k+1 endif ! Set matrix values using local numbering, defined earlier in main routine call MatSetValuesLocal(Hessian,1,row,k,col,v,INSERT_VALUES, & & info) enddo enddo ! restore vectors call VecRestoreArray(localX,x_v,x_i,info) call VecRestoreArray(Left,left_v,left_i,info) call VecRestoreArray(Right,right_v,right_i,info) call VecRestoreArray(Top,top_v,top_i,info) call VecRestoreArray(Bottom,bottom_v,bottom_i,info) ! Assemble the matrix call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,info) call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,info) call PetscLogFlops(199*xm*ym,info) return end ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h ! ---------------------------------------------------------------------------- ! !/* ! MSA_BoundaryConditions - calculates the boundary conditions for the region ! ! !*/ subroutine MSA_BoundaryConditions(info) implicit none ! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h #include "plate2f.h" integer i,j,k,limit,info,maxits integer xs, xm, gxs, gxm, ys, ym, gys, gym integer bsize, lsize, tsize, rsize PetscScalar one, two, three, tol PetscScalar scl,fnorm,det,xt,yt,hx,hy PetscScalar u1,u2,nf1,nf2,njac11,njac12,njac21,njac22 PetscScalar b, t, l, r PetscScalar boundary_v(0:1) PetscOffset boundary_i logical exitloop TaoTruth flg limit=0 maxits = 5 tol=1e-10 b=-0.5d0 t= 0.5d0 l=-0.5d0 r= 0.5d0 xt=0 yt=0 one=1.0d0 two=2.0d0 three=3.0d0 call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, & & PETSC_NULL_INTEGER,info) call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, & & gxm,gym,PETSC_NULL_INTEGER,info) bsize = xm + 2 lsize = ym + 2 rsize = ym + 2 tsize = xm + 2 call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,info) call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,info) call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,info) call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,info) hx= (r-l)/(mx+1) hy= (t-b)/(my+1) do j=0,3 if (j.eq.0) then yt=b xt=l+hx*xs limit=bsize call VecGetArray(Bottom,boundary_v,boundary_i,info) elseif (j.eq.1) then yt=t xt=l+hx*xs limit=tsize call VecGetArray(Top,boundary_v,boundary_i,info) elseif (j.eq.2) then yt=b+hy*ys xt=l limit=lsize call VecGetArray(Left,boundary_v,boundary_i,info) elseif (j.eq.3) then yt=b+hy*ys xt=r limit=rsize call VecGetArray(Right,boundary_v,boundary_i,info) endif do i=0,limit-1 u1=xt u2=-yt k = 0 exitloop = .false. do while (k .lt. maxits .and. (.not. exitloop) ) 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 .gt. tol) then 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 else exitloop = .true. endif k=k+1 enddo boundary_v(i + boundary_i) = u1*u1-u2*u2 if ((j .eq. 0) .or. (j .eq. 1)) then xt = xt + hx else yt = yt + hy endif enddo if (j.eq.0) then call VecRestoreArray(Bottom,boundary_v,boundary_i,info) elseif (j.eq.1) then call VecRestoreArray(Top,boundary_v,boundary_i,info) elseif (j.eq.2) then call VecRestoreArray(Left,boundary_v,boundary_i,info) elseif (j.eq.3) then call VecRestoreArray(Right,boundary_v,boundary_i,info) endif enddo ! Scale the boundary if desired call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom", & & scl,flg,info) if (flg .eq. PETSC_TRUE) then call VecScale(scl,Bottom,info) endif call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top", & & scl,flg,info) if (flg .eq. PETSC_TRUE) then call VecScale(scl,Top,info) endif call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right", & & scl,flg,info) if (flg .eq. PETSC_TRUE) then call VecScale(scl,Right,info) endif call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left", & & scl,flg,info) if (flg .eq. PETSC_TRUE) then call VecScale(scl,Left,info) endif return end ! ---------------------------------------------------------------------------- ! !/* ! MSA_Plate - Calculates an obstacle for surface to stretch over ! ! Output Parameter: !. xl - lower bound vector !. xu - upper bound vector ! !*/ subroutine MSA_Plate(xl,xu,info) implicit none ! mx,my,bmx,bmy,da,bheight defined in plate2f.h #include "plate2f.h" Vec xl,xu integer i,j,row,info integer xs, xm, ys, ym PetscScalar lb,ub ! PETSc's VecGetArray acts differently in Fortran than it does in C. ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info) ! will return an array of doubles referenced by x_array offset by x_index. ! i.e., to reference the kth element of X, use x_array(k + x_index). ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. PetscScalar xl_v(0:1) PetscOffset xl_i lb = -1.0d300 ub = 1.0d300 if (bmy .lt. 0) bmy = 0 if (bmy .gt. my) bmy = my if (bmx .lt. 0) bmx = 0 if (bmx .gt. mx) bmx = mx call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, & & PETSC_NULL_INTEGER,info) call VecSet(xl,lb,info) call VecSet(xu,ub,info) call VecGetArray(xl,xl_v,xl_i,info) do i=xs,xs+xm-1 do j=ys,ys+ym-1 row=(j-ys)*xm + (i-xs) if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. & & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then xl_v(xl_i+row) = bheight endif enddo enddo call VecRestoreArray(xl,xl_v,xl_i,info) return end ! ---------------------------------------------------------------------------- ! !/* ! MSA_InitialPoint - Calculates an obstacle for surface to stretch over ! ! Output Parameter: !. X - vector for initial guess ! !*/ subroutine MSA_InitialPoint(X, info) implicit none ! mx,my,localX,da,Top,Left,Bottom,Right defined in plate2f.h #include "plate2f.h" Vec X integer start,i,j,info integer row,xs,xm,gxs,gxm,ys,ym,gys,gym PetscScalar zero, np5 ! PETSc's VecGetArray acts differently in Fortran than it does in C. ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (integer) x_index, info) ! will return an array of doubles referenced by x_array offset by x_index. ! i.e., to reference the kth element of X, use x_array(k + x_index). ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. PetscScalar left_v(0:1), right_v(0:1), top_v(0:1) PetscScalar x_v(0:1), bottom_v(0:1) PetscOffset left_i, right_i, top_i, bottom_i, x_i PetscTruth flg PetscRandom rctx zero = 0.0d0 np5 = -0.5d0 call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start", & & start,flg,info) if ((flg .eq. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable call VecSet(X,zero,info) elseif ((flg .eq. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5 call PetscRandomCreate(MPI_COMM_WORLD,RANDOM_DEFAULT,rctx,info) do i=0,start-1 call VecSetRandom(X,rctx,info) enddo call PetscRandomDestroy(rctx,info) call VecShift(X,np5,info) else ! average of boundary conditions ! Get Local mesh boundaries call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, & & PETSC_NULL_INTEGER,info) call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, & & PETSC_NULL_INTEGER,info) ! Get pointers to vector data call VecGetArray(Top,top_v,top_i,info) call VecGetArray(Bottom,bottom_v,bottom_i,info) call VecGetArray(Left,left_v,left_i,info) call VecGetArray(Right,right_v,right_i,info) call VecGetArray(localX,x_v,x_i,info) ! Perform local computations do j=ys,ys+ym-1 do i=xs,xs+xm-1 row = (j-gys)*gxm + (i-gxs) x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my & & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + & & (i+1)*left_v(left_i+j-ys+1)/mx + & & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5 enddo enddo ! Restore vectors call VecRestoreArray(localX,x_v,x_i,info) call VecRestoreArray(Left,left_v,left_i,info) call VecRestoreArray(Top,top_v,top_i,info) call VecRestoreArray(Bottom,bottom_v,bottom_i,info) call VecRestoreArray(Right,right_v,right_i,info) call DALocalToGlobal(da,localX,INSERT_VALUES,X,info) endif return end