Actual source code: ex34f90.F
1: #define PETSC_AVOID_DECLARATIONS
2: #include include/finclude/petscall.h
4: !
5: ! Implements a nonlinear solver for a simple domain
6: ! consisting of 2 zero dimensional problems linked by
7: ! 2 one dimensional problems.
8: !
9: ! Channel1
10: ! -------------------------
11: ! Pool 1 Pool 2
12: ! -------------------------
13: ! Channel2
14: !
15: ! For Newton's method with block Jacobi (one block for
16: ! each channel and one block for each pool) run with
17: !
18: ! -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
19: ! -help lists all options
21: program ex34f90
22: use mex34f90
25: DMMG dmmg
26: PetscErrorCode ierr
27: DA da
28: DMComposite dm
29: type(appctx) app
30: external FormInitialGuess,FormFunction
33: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
35: ! Create the composite DM object that manages the grid
37: call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr)
38: call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1, &
39: & PETSC_NULL_INTEGER,da,ierr)
40: call DMCompositeAddDA(dm,da,ierr)
41: call DADestroy(da,ierr)
42: call DMCompositeAddArray(dm,0,app%np,ierr)
43: call DMCompositeAddDA(dm,da,ierr)
44: call DMCompositeAddArray(dm,0,app%np,ierr)
46: ! Create solver object and associate it with the unknowns (on the grid)
48: call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr)
49: call DMMGSetDM(dmmg,dm,ierr)
50: call DMCompositeDestroy(dm,ierr)
51: call DMMGSetUser(dmmg,0,app,ierr) ! currently only one level solver
52: call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr)
53: call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,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