Actual source code: ex12f.F90
1: !
2: !
3: ! This example demonstrates basic use of the SNES Fortran interface.
4: !
5: !
6: #include <petsc/finclude/petscsnes.h>
7: module ex12fmodule
8: use petscsnes
9: implicit none
10: type User
11: DM da
12: Vec F
13: Vec xl
14: MPIU_Comm comm
15: PetscInt N
16: end type User
17: type monctx
18: PetscInt :: its, lag
19: end type monctx
21: contains
22: ! -------------------- Evaluate Function F(x) ---------------------
24: subroutine FormFunction(snes, x, f, ctx, ierr)
25: SNES snes
26: Vec x, f
27: type(User) ctx
28: PetscMPIInt rank, size, zero
29: PetscInt i, s, n
30: PetscErrorCode ierr
31: PetscScalar h, d
32: PetscScalar, pointer :: vf2(:), vxx(:), vff(:)
34: zero = 0
35: PetscCallMPI(MPI_Comm_rank(ctx%comm, rank, ierr))
36: PetscCallMPI(MPI_Comm_size(ctx%comm, size, ierr))
37: h = 1.0/(real(ctx%N) - 1.0)
38: PetscCall(DMGlobalToLocalBegin(ctx%da, x, INSERT_VALUES, ctx%xl, ierr))
39: PetscCall(DMGlobalToLocalEnd(ctx%da, x, INSERT_VALUES, ctx%xl, ierr))
41: PetscCall(VecGetLocalSize(ctx%xl, n, ierr))
42: if (n > 1000) then
43: print *, 'Local work array not big enough'
44: call MPI_Abort(PETSC_COMM_WORLD, zero, ierr)
45: end if
47: PetscCall(VecGetArrayRead(ctx%xl, vxx, ierr))
48: PetscCall(VecGetArray(f, vff, ierr))
49: PetscCall(VecGetArray(ctx%F, vF2, ierr))
51: d = h*h
53: !
54: ! Note that the array vxx() was obtained from a ghosted local vector
55: ! ctx%xl while the array vff() was obtained from the non-ghosted parallel
56: ! vector F. This is why there is a need for shift variable s. Since vff()
57: ! does not have locations for the ghost variables we need to index in it
58: ! slightly different then indexing into vxx(). For example on processor
59: ! 1 (the second processor)
60: !
61: ! xx(1) xx(2) xx(3) .....
62: ! ^^^^^^^ ^^^^^ ^^^^^
63: ! ghost value 1st local value 2nd local value
64: !
65: ! ff(1) ff(2)
66: ! ^^^^^^^ ^^^^^^^
67: ! 1st local value 2nd local value
68: !
69: if (rank == 0) then
70: s = 0
71: vff(1) = vxx(1)
72: else
73: s = 1
74: end if
76: do i = 1, n - 2
77: vff(i - s + 1) = d*(vxx(i) - 2.0*vxx(i + 1) + vxx(i + 2)) + vxx(i + 1)*vxx(i + 1) - vF2(i - s + 1)
78: end do
80: if (rank == size - 1) then
81: vff(n - s) = vxx(n) - 1.0
82: end if
84: PetscCall(VecRestoreArray(f, vff, ierr))
85: PetscCall(VecRestoreArrayRead(ctx%xl, vxx, ierr))
86: PetscCall(VecRestoreArray(ctx%F, vF2, ierr))
87: end
89: ! ---------------------------------------------------------------------
90: ! Subroutine FormMonitor
91: ! This function lets up keep track of the SNES progress at each step
92: ! In this routine, we determine when the Jacobian is rebuilt with the parameter 'jag'
93: !
94: ! Input Parameters:
95: ! snes - SNES nonlinear solver context
96: ! its - current nonlinear iteration, starting from a call of SNESSolve()
97: ! norm - 2-norm of current residual (may be approximate)
98: ! snesm - monctx designed module (included in Snesmmod)
99: ! ---------------------------------------------------------------------
100: subroutine FormMonitor(snes, its, norm, snesm, ierr)
102: SNES :: snes
103: PetscInt :: its, one, mone
104: PetscScalar :: norm
105: type(monctx) :: snesm
106: PetscErrorCode :: ierr
108: ! write(6,*) ' '
109: ! write(6,*) ' its ',its,snesm%its,'lag',
110: ! & snesm%lag
111: ! call flush(6)
112: if (mod(snesm%its, snesm%lag) == 0) then
113: one = 1
114: PetscCall(SNESSetLagJacobian(snes, one, ierr)) ! build jacobian
115: else
116: mone = -1
117: PetscCall(SNESSetLagJacobian(snes, mone, ierr)) ! do NOT build jacobian
118: end if
119: snesm%its = snesm%its + 1
120: end subroutine FormMonitor
122: ! -------------------- Form initial approximation -----------------
124: subroutine FormInitialGuess(snes, x, ierr)
126: PetscErrorCode ierr
127: Vec x
128: SNES snes
129: PetscScalar five
131: five = .5
132: PetscCall(VecSet(x, five, ierr))
133: end
135: ! -------------------- Evaluate Jacobian --------------------
137: subroutine FormJacobian(snes, x, jac, B, ctx, ierr)
139: SNES snes
140: Vec x
141: Mat jac, B
142: type(User) ctx
143: PetscInt ii, istart, iend
144: PetscInt i, j, n, end, start, i1
145: PetscErrorCode ierr
146: PetscMPIInt rank, size
147: PetscScalar d, A, h
148: PetscScalar, pointer :: vxx(:)
150: i1 = 1
151: h = 1.0/(real(ctx%N) - 1.0)
152: d = h*h
153: PetscCallMPI(MPI_Comm_rank(ctx%comm, rank, ierr))
154: PetscCallMPI(MPI_Comm_size(ctx%comm, size, ierr))
156: PetscCall(VecGetArrayRead(x, vxx, ierr))
157: PetscCall(VecGetOwnershipRange(x, start, end, ierr))
158: n = end - start
160: if (rank == 0) then
161: A = 1.0
162: PetscCall(MatSetValues(jac, i1, [start], i1, [start], [A], INSERT_VALUES, ierr))
163: istart = 1
164: else
165: istart = 0
166: end if
167: if (rank == size - 1) then
168: i = INT(ctx%N - 1)
169: A = 1.0
170: PetscCall(MatSetValues(jac, i1, [i], i1, [i], [A], INSERT_VALUES, ierr))
171: iend = n - 1
172: else
173: iend = n
174: end if
175: do i = istart, iend - 1
176: ii = i + start
177: j = start + i - 1
178: PetscCall(MatSetValues(jac, i1, [ii], i1, [j], [d], INSERT_VALUES, ierr))
179: j = start + i + 1
180: PetscCall(MatSetValues(jac, i1, [ii], i1, [j], [d], INSERT_VALUES, ierr))
181: A = -2.0*d + 2.0*vxx(i + 1)
182: PetscCall(MatSetValues(jac, i1, [ii], i1, [ii], [A], INSERT_VALUES, ierr))
183: end do
184: PetscCall(VecRestoreArrayRead(x, vxx, ierr))
185: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
186: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
187: end
189: end module
191: program main
192: use ex12fmodule
193: implicit none
194: type(User) ctx
195: PetscMPIInt rank, size
196: PetscErrorCode ierr
197: PetscInt N, start, end, nn, i
198: PetscInt ii, its, i1, i0, i3
199: PetscBool flg
200: SNES snes
201: Mat J
202: Vec x, r, u
203: PetscScalar xp, FF, UU, h
204: character*(10) matrixname
205: type(monctx) :: snesm
207: PetscCallA(PetscInitialize(ierr))
208: i1 = 1
209: i0 = 0
210: i3 = 3
211: N = 10
212: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', N, flg, ierr))
213: h = 1.0/real(N - 1)
214: ctx%N = N
215: ctx%comm = PETSC_COMM_WORLD
217: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
218: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
220: ! Set up data structures
221: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, i1, i1, PETSC_NULL_INTEGER_ARRAY, ctx%da, ierr))
222: PetscCallA(DMSetFromOptions(ctx%da, ierr))
223: PetscCallA(DMSetUp(ctx%da, ierr))
224: PetscCallA(DMCreateGlobalVector(ctx%da, x, ierr))
225: PetscCallA(DMCreateLocalVector(ctx%da, ctx%xl, ierr))
227: PetscCallA(PetscObjectSetName(x, 'Approximate Solution', ierr))
228: PetscCallA(VecDuplicate(x, r, ierr))
229: PetscCallA(VecDuplicate(x, ctx%F, ierr))
230: PetscCallA(VecDuplicate(x, U, ierr))
231: PetscCallA(PetscObjectSetName(U, 'Exact Solution', ierr))
233: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, i3, PETSC_NULL_INTEGER_ARRAY, i0, PETSC_NULL_INTEGER_ARRAY, J, ierr))
234: PetscCallA(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE, ierr))
235: PetscCallA(MatGetType(J, matrixname, ierr))
237: ! Store right-hand side of PDE and exact solution
238: PetscCallA(VecGetOwnershipRange(x, start, end, ierr))
239: xp = h*start
240: nn = end - start
241: ii = start
242: do i = 0, nn - 1
243: FF = 6.0*xp + (xp + 1.e-12)**6.e0
244: UU = xp*xp*xp
245: PetscCallA(VecSetValues(ctx%F, i1, [ii], [FF], INSERT_VALUES, ierr))
246: PetscCallA(VecSetValues(U, i1, [ii], [UU], INSERT_VALUES, ierr))
247: xp = xp + h
248: ii = ii + 1
249: end do
250: PetscCallA(VecAssemblyBegin(ctx%F, ierr))
251: PetscCallA(VecAssemblyEnd(ctx%F, ierr))
252: PetscCallA(VecAssemblyBegin(U, ierr))
253: PetscCallA(VecAssemblyEnd(U, ierr))
255: ! Create nonlinear solver
256: PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
258: ! Set various routines and options
259: PetscCallA(SNESSetFunction(snes, r, FormFunction, ctx, ierr))
260: PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, ctx, ierr))
262: snesm%its = 0
263: PetscCallA(SNESGetLagJacobian(snes, snesm%lag, ierr))
264: PetscCallA(SNESMonitorSet(snes, FormMonitor, snesm, PETSC_NULL_FUNCTION, ierr))
265: PetscCallA(SNESSetFromOptions(snes, ierr))
267: ! Solve nonlinear system
268: PetscCallA(FormInitialGuess(snes, x, ierr))
269: PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
270: PetscCallA(SNESGetIterationNumber(snes, its, ierr))
272: ! Free work space. All PETSc objects should be destroyed when they
273: ! are no longer needed.
274: PetscCallA(VecDestroy(x, ierr))
275: PetscCallA(VecDestroy(ctx%xl, ierr))
276: PetscCallA(VecDestroy(r, ierr))
277: PetscCallA(VecDestroy(U, ierr))
278: PetscCallA(VecDestroy(ctx%F, ierr))
279: PetscCallA(MatDestroy(J, ierr))
280: PetscCallA(SNESDestroy(snes, ierr))
281: PetscCallA(DMDestroy(ctx%da, ierr))
282: PetscCallA(PetscFinalize(ierr))
283: end
285: !/*TEST
286: !
287: ! test:
288: ! nsize: 2
289: ! args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
290: ! output_file: output/ex12_1.out
291: !
292: !TEST*/