Actual source code: ex34f90.F

  1: #include "finclude/petscdef.h"

  3: !
  4: !      Implements a nonlinear solver for a simple domain
  5: !     consisting of 2 zero dimensional problems linked by
  6: !     2 one dimensional problems.
  7: !
  8: !                           Channel1
  9: !                       -------------------------
 10: !               Pool 1                              Pool 2
 11: !                       -------------------------
 12: !                           Channel2
 13: !
 14: !     For Newton's method with block Jacobi (one block for
 15: !     each channel and one block for each pool) run with
 16: !
 17: !       -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
 18: !       -help lists all options

 20:       program ex34f90
 21:       use mex34f90


 24:       DMMG             dmmg
 25:       PetscErrorCode   ierr
 26:       DA               da
 27:       DMComposite      dm
 28:       type(appctx)     app
 29:       external         FormInitialGuess,FormFunction


 32:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 34: !      Create the composite DM object that manages the grid

 36:       call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr)
 37:       call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1, &
 38:      &                PETSC_NULL_INTEGER,da,ierr)
 39:       call DMCompositeAddDM(dm,da,ierr)
 40:       call DADestroy(da,ierr)
 41:       call DMCompositeAddArray(dm,0,app%np,ierr)
 42:       call DMCompositeAddDM(dm,da,ierr)
 43:       call DMCompositeAddArray(dm,0,app%np,ierr)

 45: !       Create solver object and associate it with the unknowns (on the grid)

 47:       call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr)
 48:       call DMMGSetDM(dmmg,dm,ierr)
 49:       call DMCompositeDestroy(dm,ierr)
 50:       call DMMGSetUser(dmmg,0,app,ierr)  ! currently only one level solver
 51:       call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr)
 52:       call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,ierr)
 53:       call DMMGSetFromOptions(dmmg,ierr)

 55: !      Solve the nonlinear system
 56: !
 57:       call DMMGSolve(dmmg,ierr)

 59:       call DMMGDestroy(dmmg,ierr)
 60:       call PetscFinalize(ierr)
 61:       end


 64:       subroutine FormInitialGuess(dmmg,v,ierr)
 65:       use mex34f90
 66:       use mex34finterfaces

 68:       DMMG dmmg
 69:       Vec v
 70:       PetscErrorCode ierr

 72:       DMComposite dm
 73:       Vec  vc1,vc2
 74:       PetscInt np,i
 75:       DA da
 76:       type(DALocalInfof90) dainfo
 77:       type(poolfield), pointer :: Pool1,Pool2
 78:       type(channelfield), pointer :: v1(:),v2(:)
 79:       type(appctx), pointer :: app
 80:       PetscMPIInt rank

 82:       call DMMGGetUser(dmmg,app,ierr)   ! access user context

 84: !       Access the grid information

 86:       call DMMGGetDM(dmmg,dm,ierr)
 87:       call DMCompositeGetEntries(dm,da,np,da,np,ierr) ! get the da's and sizes that define the unknowns
 88:       call DAGetLocalInfof90(da,dainfo,ierr)

 90: !      Acess the vector unknowns in global (nonghosted) form

 92:       call DMCompositeGetAccess(dm,v,vc1,Pool1,vc2,Pool2,ierr)  !
 93:       call DAVecGetArrayf90(da,vc1,v1,ierr)
 94:       call DAVecGetArrayf90(da,vc2,v2,ierr)

 96:       call MPI_Comm_rank(app%comm,rank,ierr)
 97: !      the pools are only unknowns on process 0
 98:       if (rank == 0) then
 99:          Pool1%p0 = -2.0
100:          Pool1%p1 = -2000.0
101:          Pool2%p0 = -20
102:          Pool2%p1 = -200
103:       endif

105: !       put some random numbers as an initial guess
106:       do i=dainfo%xs,dainfo%xs+dainfo%xm-1
107:         v1(i)%c0 = 3*i - .5
108:         v1(i)%c1 = 3*i - (1.d0/3.d0)
109:         v1(i)%c2 = 3*i - .1
110:       enddo

112:       call DAVecRestoreArrayf90(da,vc1,v1,ierr)
113:       call DAVecRestoreArrayf90(da,vc2,v2,ierr)
114:       call DMCompositeRestoreAccess(dm,v,vc1,Pool1,vc2,Pool2,ierr)

116:       CHKMEMQ
117:       return
118:       end

120:       subroutine FormFunction(snes,x,f,dmmg,ierr)
121:       use mex34f90
122:       use mex34finterfaces

124:       SNES snes
125:       Vec x,f
126:       DMMG dmmg
127:       PetscErrorCode ierr

129:       DMComposite dm
130:       Vec  fvc1,fvc2,xvc1,xvc2
131:       PetscInt np,i
132:       DA da
133:       type(DALocalInfof90) dainfo
134:       type(poolfield), pointer :: fPool1,fPool2
135:       type(poolfield), pointer :: xPool1,xPool2
136:       type(channelfield), pointer :: fv1(:),fv2(:),xv1(:),xv2(:)
137:       type(appctx), pointer :: app
138:       PetscMPIInt rank

140:       call DMMGGetUser(dmmg,app,ierr)   ! access user context

142: !         Access the grid information

144:       call DMMGGetDM(dmmg,dm,ierr)
145:       call DMCompositeGetEntries(dm,da,np,da,np,ierr) ! get the da's and sizes that define the unknowns
146:       call DAGetLocalInfof90(da,dainfo,ierr)

148: !        Access the local (ghosted) parts of x

150:       call DMCompositeGetLocalVectors(dm,xvc1,xPool1,xvc2,xPool2,ierr)
151:       call DMCompositeScatter(dm,x,xvc1,xPool1,xvc2,xPool2,ierr)
152:       call DAVecGetArrayf90(da,xvc1,xv1,ierr)
153:       call DAVecGetArrayf90(da,xvc2,xv2,ierr)

155: !       Access the global (non-ghosted) parts of f

157:       call DMCompositeGetAccess(dm,f,fvc1,fPool1,fvc2,fPool2,ierr)
158:       call DAVecGetArrayf90(da,fvc1,fv1,ierr)
159:       call DAVecGetArrayf90(da,fvc2,fv2,ierr)


162:       call MPI_Comm_rank(app%comm,rank,ierr)
163: !      fPool only lives on zeroth process, ghosted xPool lives on all processes
164:       if (rank == 0) then
165:          fPool1%p0 = xPool1%p0
166:          fPool1%p1 = xPool1%p1
167:          fPool2%p0 = xPool2%p0
168:          fPool2%p1 = xPool2%p1
169:       endif

171: !        Replace with some difference scheme
172:       do i=dainfo%xs,dainfo%xs+dainfo%xm-1
173:         fv1(i)%c0 = xv1(i)%c0
174:         fv1(i)%c1 = xv1(i)%c1
175:         fv1(i)%c2 = xv1(i)%c2

177:         fv2(i)%c0 = xv2(i)%c0
178:         fv2(i)%c1 = xv2(i)%c1
179:         fv2(i)%c2 = xv2(i)%c2
180:       enddo

182:       call DAVecRestoreArrayf90(da,xvc1,xv1,ierr)
183:       call DAVecRestoreArrayf90(da,xvc2,xv2,ierr)
184:       call DMCompositeRestoreLocalVectors(dm,xvc1,xPool1,xvc2,xPool2,   &
185:      &     ierr)

187:       call DAVecRestoreArrayf90(da,fvc1,fv1,ierr)
188:       call DAVecRestoreArrayf90(da,fvc2,fv2,ierr)
189:       call DMCompositeRestoreAccess(dm,f,fvc1,fPool1,fvc2,fPool2,ierr)

191:       return
192:       end