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