Actual source code: ex1f.F90
1: !
2: ! Solves the time dependent Bratu problem using pseudo-timestepping
3: !
4: ! This code demonstrates how one may solve a nonlinear problem
5: ! with pseudo-timestepping. In this simple example, the pseudo-timestep
6: ! is the same for all grid points, i.e., this is equivalent to using
7: ! the backward Euler method with a variable timestep.
8: !
9: ! Note: This example does not require pseudo-timestepping since it
10: ! is an easy nonlinear problem, but it is included to demonstrate how
11: ! the pseudo-timestepping may be done.
12: !
13: ! See snes/tutorials/ex4.c[ex4f.F] and
14: ! snes/tutorials/ex5.c[ex5f.F] where the problem is described
15: ! and solved using the method of Newton alone.
16: !
17: !
18: !23456789012345678901234567890123456789012345678901234567890123456789012
19: #include <petsc/finclude/petscts.h>
20: module ex1fmodule
21: use petscts
22: implicit none
23: PetscInt, parameter :: param = 1, lmx = 2, lmy = 3
24: contains
25: !
26: ! -------------------- Evaluate Function F(x) ---------------------
27: !
28: subroutine FormFunction(ts, t, X, F, user, ierr)
30: TS ts
31: PetscReal t
32: Vec X, F
33: PetscReal user(3)
34: PetscErrorCode, intent(out) :: ierr
35: PetscInt i, j, row, mx, my
36: PetscReal lambda
37: PetscReal hx, hy, hxdhy, hydhx
38: PetscScalar ut, ub, ul, ur, u
39: PetscScalar uxx, uyy, sc
40: PetscScalar, pointer :: xx(:), ff(:)
42: mx = int(user(lmx))
43: my = int(user(lmy))
44: lambda = user(param)
46: hx = 1.0/real(mx - 1)
47: hy = 1.0/real(my - 1)
48: sc = hx*hy
49: hxdhy = hx/hy
50: hydhx = hy/hx
52: PetscCall(VecGetArrayRead(X, xx, ierr))
53: PetscCall(VecGetArray(F, ff, ierr))
54: do j = 1, my
55: do i = 1, mx
56: row = i + (j - 1)*mx
57: if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
58: ff(row) = xx(row)
59: else
60: u = xx(row)
61: ub = xx(row - mx)
62: ul = xx(row - 1)
63: ut = xx(row + mx)
64: ur = xx(row + 1)
65: uxx = (-ur + 2.0_PETSC_REAL_KIND*u - ul)*hydhx
66: uyy = (-ut + 2.0_PETSC_REAL_KIND*u - ub)*hxdhy
67: ff(row) = -uxx - uyy + sc*lambda*exp(u)
68: end if
69: end do
70: end do
72: PetscCall(VecRestoreArrayRead(X, xx, ierr))
73: PetscCall(VecRestoreArray(F, ff, ierr))
74: end
75: !
76: ! -------------------- Evaluate Jacobian of F(x) --------------------
77: !
78: subroutine FormJacobian(ts, ctime, X, JJ, B, user, ierr)
80: TS ts
81: Vec X
82: Mat JJ, B
83: PetscReal user(3), ctime
84: PetscErrorCode, intent(out) :: ierr
85: Mat jac
86: PetscInt i, j, row(1), mx, my
87: PetscInt col(5)
88: PetscScalar lambda
89: PetscScalar v(5), sc
90: PetscScalar, pointer :: xx(:)
91: PetscReal hx, hy, hxdhy, hydhx
93: jac = B
95: mx = int(user(lmx))
96: my = int(user(lmy))
97: lambda = user(param)
99: hx = 1.0/real(mx - 1)
100: hy = 1.0/real(my - 1)
101: sc = hx*hy
102: hxdhy = hx/hy
103: hydhx = hy/hx
105: PetscCall(VecGetArrayRead(X, xx, ierr))
106: do j = 1, my
107: do i = 1, mx
108: !
109: ! When inserting into PETSc matrices, indices start at 0
110: !
111: row(1) = i - 1 + (j - 1)*mx
112: if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
113: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [row], 1_PETSC_INT_KIND, [row], [1.0_PETSC_REAL_KIND], INSERT_VALUES, ierr))
114: else
115: v(1) = hxdhy
116: col(1) = row(1) - mx
117: v(2) = hydhx
118: col(2) = row(1) - 1
119: v(3) = -2.0_PETSC_REAL_KIND*(hydhx + hxdhy) + sc*lambda*exp(xx(row(1)))
120: col(3) = row(1)
121: v(4) = hydhx
122: col(4) = row(1) + 1
123: v(5) = hxdhy
124: col(5) = row(1) + mx
125: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [row], 5_PETSC_INT_KIND, col, v, INSERT_VALUES, ierr))
126: end if
127: end do
128: end do
129: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
130: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
131: PetscCall(VecRestoreArray(X, xx, ierr))
132: end
133: !
134: ! -------------------- Form initial approximation -----------------
135: !
136: subroutine FormInitialGuess(X, user, ierr)
138: Vec X
139: PetscReal user(3)
140: PetscInt i, j, row, mx, my
141: PetscErrorCode, intent(out) :: ierr
142: PetscReal lambda
143: PetscReal temp1, temp, hx, hy
144: PetscScalar, pointer :: xx(:)
146: mx = int(user(lmx))
147: my = int(user(lmy))
148: lambda = user(param)
150: hy = 1._PETSC_REAL_KIND/(my - 1)
151: hx = 1._PETSC_REAL_KIND/(mx - 1)
153: PetscCall(VecGetArray(X, xx, ierr))
154: temp1 = lambda/(lambda + 1._PETSC_REAL_KIND)
155: do j = 1, my
156: temp = min(j - 1, my - j)*hy
157: do i = 1, mx
158: row = i + (j - 1)*mx
159: if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then
160: xx(row) = 0.0
161: else
162: xx(row) = temp1*sqrt(min(min(i - 1, mx - i)*hx, temp))
163: end if
164: end do
165: end do
166: PetscCall(VecRestoreArray(X, xx, ierr))
167: end
168: end module ex1fmodule
170: program main
171: use petscts
172: use ex1fmodule
173: implicit none
174: !
175: ! Create an application context to contain data needed by the
176: ! application-provided call-back routines, FormJacobian() and
177: ! FormFunction(). We use a double precision array with three
178: ! entries indexed by param, lmx, lmy.
179: !
180: PetscReal user(3)
181: !
182: ! Data for problem
183: !
184: TS ts
185: Vec x, r
186: Mat J
187: PetscInt its, N, itmp
188: PetscBool flg
189: PetscErrorCode ierr
190: PetscReal, parameter :: param_min = 0.0, param_max = 6.81, tmax = 1.e12
191: PetscReal dt, ftime
192: TSConvergedReason reason
194: PetscCallA(PetscInitialize(ierr))
195: user(lmx) = 4.0
196: user(lmy) = 4.0
197: user(param) = 6.0
199: !
200: ! Allow user to set the grid dimensions and nonlinearity parameter at run-time
201: !
202: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', user(lmx), flg, ierr))
203: itmp = 4
204: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', itmp, flg, ierr))
205: user(lmy) = itmp
206: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-param', user(param), flg, ierr))
207: if (user(param) >= param_max .or. user(param) <= param_min) then
208: print *, 'Parameter is out of range'
209: end if
210: if (user(lmx) > user(lmy)) then
211: dt = .5/user(lmx)
212: else
213: dt = .5/user(lmy)
214: end if
215: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-dt', dt, flg, ierr))
216: N = int(user(lmx)*user(lmy))
218: ! Create vectors to hold the solution and function value
219: PetscCallA(VecCreateSeq(PETSC_COMM_SELF, N, x, ierr))
220: PetscCallA(VecDuplicate(x, r, ierr))
222: ! Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
223: ! in the sparse matrix. Note that this is not the optimal strategy see
224: ! the Performance chapter of the users manual for information on
225: ! preallocating memory in sparse matrices.
226: PetscCallA(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 5_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, J, ierr))
228: ! Create timestepper context
229: PetscCallA(TSCreate(PETSC_COMM_WORLD, ts, ierr))
230: PetscCallA(TSSetProblemType(ts, TS_NONLINEAR, ierr))
232: ! Tell the timestepper context where to compute solutions
233: PetscCallA(TSSetSolution(ts, x, ierr))
235: ! Provide the call-back for the nonlinear function we are
236: ! evaluating. Thus whenever the timestepping routines need the
237: ! function they will call this routine. Note the final argument
238: ! is the application context used by the call-back functions.
239: PetscCallA(TSSetRHSFunction(ts, PETSC_NULL_VEC, FormFunction, user, ierr))
241: ! Set the Jacobian matrix and the function used to compute
242: ! Jacobians.
243: PetscCallA(TSSetRHSJacobian(ts, J, J, FormJacobian, user, ierr))
245: ! For the initial guess for the problem
246: PetscCallA(FormInitialGuess(x, user, ierr))
248: ! This indicates that we are using pseudo timestepping to
249: ! find a steady state solution to the nonlinear problem.
250: PetscCallA(TSSetType(ts, TSPSEUDO, ierr))
252: ! Set the initial time to start at (this is arbitrary for
253: ! steady state problems and the initial timestep given above
254: PetscCallA(TSSetTimeStep(ts, dt, ierr))
256: ! Set a large number of timesteps and final duration time
257: ! to insure convergence to steady state.
258: PetscCallA(TSSetMaxSteps(ts, 1000_PETSC_INT_KIND, ierr))
259: PetscCallA(TSSetMaxTime(ts, tmax, ierr))
260: PetscCallA(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER, ierr))
262: ! Set any additional options from the options database. This
263: ! includes all options for the nonlinear and linear solvers used
264: ! internally the timestepping routines.
265: PetscCallA(TSSetFromOptions(ts, ierr))
267: PetscCallA(TSSetUp(ts, ierr))
269: ! Perform the solve. This is where the timestepping takes place.
270: PetscCallA(TSSolve(ts, x, ierr))
271: PetscCallA(TSGetSolveTime(ts, ftime, ierr))
272: PetscCallA(TSGetStepNumber(ts, its, ierr))
273: PetscCallA(TSGetConvergedReason(ts, reason, ierr))
275: write (6, 100) its, ftime, reason
276: 100 format('Number of pseudo time-steps ', i5, ' final time ', 1pe9.2, ' reason ', i3)
278: !
279: ! Free the data structures constructed above
280: !
281: PetscCallA(VecDestroy(x, ierr))
282: PetscCallA(VecDestroy(r, ierr))
283: PetscCallA(MatDestroy(J, ierr))
284: PetscCallA(TSDestroy(ts, ierr))
285: PetscCallA(PetscFinalize(ierr))
286: end
287: !/*TEST
288: !
289: ! test:
290: ! TODO: broken
291: ! args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_monitor_pseudo -ts_max_snes_failures 3 -ts_pseudo_frtol 1.e-5 -snes_stol 1e-5
292: !
293: !TEST*/