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, parameter :: zero = 0
29: PetscMPIInt rank, size
30: PetscInt i, s, n
31: PetscErrorCode ierr
32: PetscScalar h, d
33: PetscScalar, pointer :: vf2(:), vxx(:), vff(:)
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**2
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
145: PetscErrorCode, intent(out) :: ierr
146: PetscMPIInt rank, size
147: PetscScalar d, A, h
148: PetscScalar, pointer :: vxx(:)
150: h = 1.0/(real(ctx%N) - 1.0)
151: d = h**2
152: PetscCallMPI(MPI_Comm_rank(ctx%comm, rank, ierr))
153: PetscCallMPI(MPI_Comm_size(ctx%comm, size, ierr))
155: PetscCall(VecGetArrayRead(x, vxx, ierr))
156: PetscCall(VecGetOwnershipRange(x, start, end, ierr))
157: n = end - start
159: if (rank == 0) then
160: A = 1.0
161: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [start], 1_PETSC_INT_KIND, [start], [A], INSERT_VALUES, ierr))
162: istart = 1
163: else
164: istart = 0
165: end if
166: if (rank == size - 1) then
167: i = INT(ctx%N - 1)
168: A = 1.0
169: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [i], 1_PETSC_INT_KIND, [i], [A], INSERT_VALUES, ierr))
170: iend = n - 1
171: else
172: iend = n
173: end if
174: do i = istart, iend - 1
175: ii = i + start
176: j = start + i - 1
177: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [ii], 1_PETSC_INT_KIND, [j], [d], INSERT_VALUES, ierr))
178: j = start + i + 1
179: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [ii], 1_PETSC_INT_KIND, [j], [d], INSERT_VALUES, ierr))
180: A = -2.0*d + 2.0*vxx(i + 1)
181: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [ii], 1_PETSC_INT_KIND, [ii], [A], INSERT_VALUES, ierr))
182: end do
183: PetscCall(VecRestoreArrayRead(x, vxx, ierr))
184: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
185: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
186: end
188: end module
190: program main
191: use ex12fmodule
192: implicit none
193: type(User) ctx
194: PetscMPIInt rank, size
195: PetscErrorCode ierr
196: PetscInt N, start, end, nn, i
197: PetscInt ii, its
198: PetscBool flg
199: SNES snes
200: Mat J
201: Vec x, r, u
202: PetscScalar xp, FF, UU, h
203: character(len=10) matrixname
204: type(monctx) :: snesm
206: PetscCallA(PetscInitialize(ierr))
207: N = 10
208: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', N, flg, ierr))
209: h = 1.0/real(N - 1)
210: ctx%N = N
211: ctx%comm = PETSC_COMM_WORLD
213: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
214: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
216: ! Set up data structures
217: PetscCallA(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 1_PETSC_INT_KIND, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ctx%da, ierr))
218: PetscCallA(DMSetFromOptions(ctx%da, ierr))
219: PetscCallA(DMSetUp(ctx%da, ierr))
220: PetscCallA(DMCreateGlobalVector(ctx%da, x, ierr))
221: PetscCallA(DMCreateLocalVector(ctx%da, ctx%xl, ierr))
223: PetscCallA(PetscObjectSetName(x, 'Approximate Solution', ierr))
224: PetscCallA(VecDuplicate(x, r, ierr))
225: PetscCallA(VecDuplicate(x, ctx%F, ierr))
226: PetscCallA(VecDuplicate(x, U, ierr))
227: PetscCallA(PetscObjectSetName(U, 'Exact Solution', ierr))
229: PetscCallA(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 0_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, J, ierr))
230: PetscCallA(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE, ierr))
231: PetscCallA(MatGetType(J, matrixname, ierr))
233: ! Store right-hand side of PDE and exact solution
234: PetscCallA(VecGetOwnershipRange(x, start, end, ierr))
235: xp = h*start
236: nn = end - start
237: ii = start
238: do i = 0, nn - 1
239: FF = 6.0*xp + (xp + 1.e-12)**6.e0
240: UU = xp*xp*xp
241: PetscCallA(VecSetValues(ctx%F, 1_PETSC_INT_KIND, [ii], [FF], INSERT_VALUES, ierr))
242: PetscCallA(VecSetValues(U, 1_PETSC_INT_KIND, [ii], [UU], INSERT_VALUES, ierr))
243: xp = xp + h
244: ii = ii + 1
245: end do
246: PetscCallA(VecAssemblyBegin(ctx%F, ierr))
247: PetscCallA(VecAssemblyEnd(ctx%F, ierr))
248: PetscCallA(VecAssemblyBegin(U, ierr))
249: PetscCallA(VecAssemblyEnd(U, ierr))
251: ! Create nonlinear solver
252: PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr))
254: ! Set various routines and options
255: PetscCallA(SNESSetFunction(snes, r, FormFunction, ctx, ierr))
256: PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, ctx, ierr))
258: snesm%its = 0
259: PetscCallA(SNESGetLagJacobian(snes, snesm%lag, ierr))
260: PetscCallA(SNESMonitorSet(snes, FormMonitor, snesm, PETSC_NULL_FUNCTION, ierr))
261: PetscCallA(SNESSetFromOptions(snes, ierr))
263: ! Solve nonlinear system
264: PetscCallA(FormInitialGuess(snes, x, ierr))
265: PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr))
266: PetscCallA(SNESGetIterationNumber(snes, its, ierr))
268: ! Free work space. All PETSc objects should be destroyed when they
269: ! are no longer needed.
270: PetscCallA(VecDestroy(x, ierr))
271: PetscCallA(VecDestroy(ctx%xl, ierr))
272: PetscCallA(VecDestroy(r, ierr))
273: PetscCallA(VecDestroy(U, ierr))
274: PetscCallA(VecDestroy(ctx%F, ierr))
275: PetscCallA(MatDestroy(J, ierr))
276: PetscCallA(SNESDestroy(snes, ierr))
277: PetscCallA(DMDestroy(ctx%da, ierr))
278: PetscCallA(PetscFinalize(ierr))
279: end
281: !/*TEST
282: !
283: ! test:
284: ! nsize: 2
285: ! args: -ksp_gmres_cgs_refinement_type refine_always -n 10 -snes_monitor_short
286: ! output_file: output/ex12_1.out
287: !
288: !TEST*/