Actual source code: ex1f.F
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: program main
20: #include <petsc/finclude/petscts.h>
21: use petscts
22: implicit none
24: !
25: ! Create an application context to contain data needed by the
26: ! application-provided call-back routines, FormJacobian() and
27: ! FormFunction(). We use a double precision array with three
28: ! entries indexed by param, lmx, lmy.
29: !
30: PetscReal user(3)
31: PetscInt param,lmx,lmy,i5
32: parameter (param = 1,lmx = 2,lmy = 3)
33: !
34: ! User-defined routines
35: !
36: external FormJacobian,FormFunction
37: !
38: ! Data for problem
39: !
40: TS ts
41: Vec x,r
42: Mat J
43: PetscInt its,N,i1000,itmp
44: PetscBool flg
45: PetscErrorCode ierr
46: PetscReal param_max,param_min,dt
47: PetscReal tmax
48: PetscReal ftime
49: TSConvergedReason reason
51: i5 = 5
52: param_max = 6.81
53: param_min = 0
55: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
56: if (ierr .ne. 0) then
57: print*,'Unable to initialize PETSc'
58: stop
59: endif
60: user(lmx) = 4
61: user(lmy) = 4
62: user(param) = 6.0
64: !
65: ! Allow user to set the grid dimensions and nonlinearity parameter at run-time
66: !
67: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
68: & '-mx',user(lmx),flg,ierr)
69: itmp = 4
70: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
71: & '-my',itmp,flg,ierr)
72: user(lmy) = itmp
73: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
74: & '-param',user(param),flg,ierr)
75: if (user(param) .ge. param_max .or. &
76: & user(param) .le. param_min) then
77: print*,'Parameter is out of range'
78: endif
79: if (user(lmx) .gt. user(lmy)) then
80: dt = .5/user(lmx)
81: else
82: dt = .5/user(lmy)
83: endif
84: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
85: & '-dt',dt,flg,ierr)
86: N = int(user(lmx)*user(lmy))
88: !
89: ! Create vectors to hold the solution and function value
90: !
91: call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr)
92: call VecDuplicate(x,r,ierr)
94: !
95: ! Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
96: ! in the sparse matrix. Note that this is not the optimal strategy see
97: ! the Performance chapter of the users manual for information on
98: ! preallocating memory in sparse matrices.
99: !
100: i5 = 5
101: call MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,i5,PETSC_NULL_INTEGER, &
102: & J,ierr)
104: !
105: ! Create timestepper context
106: !
108: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
109: call TSSetProblemType(ts,TS_NONLINEAR,ierr)
111: !
112: ! Tell the timestepper context where to compute solutions
113: !
115: call TSSetSolution(ts,x,ierr)
117: !
118: ! Provide the call-back for the nonlinear function we are
119: ! evaluating. Thus whenever the timestepping routines need the
120: ! function they will call this routine. Note the final argument
121: ! is the application context used by the call-back functions.
122: !
124: call TSSetRHSFunction(ts,PETSC_NULL_VEC,FormFunction,user,ierr)
126: !
127: ! Set the Jacobian matrix and the function used to compute
128: ! Jacobians.
129: !
131: call TSSetRHSJacobian(ts,J,J,FormJacobian,user,ierr)
133: !
134: ! For the initial guess for the problem
135: !
137: call FormInitialGuess(x,user,ierr)
139: !
140: ! This indicates that we are using pseudo timestepping to
141: ! find a steady state solution to the nonlinear problem.
142: !
144: call TSSetType(ts,TSPSEUDO,ierr)
146: !
147: ! Set the initial time to start at (this is arbitrary for
148: ! steady state problems and the initial timestep given above
149: !
151: call TSSetTimeStep(ts,dt,ierr)
153: !
154: ! Set a large number of timesteps and final duration time
155: ! to insure convergence to steady state.
156: !
157: i1000 = 1000
158: tmax = 1.e12
159: call TSSetMaxSteps(ts,i1000,ierr)
160: call TSSetMaxTime(ts,tmax,ierr)
161: call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)
163: !
164: ! Set any additional options from the options database. This
165: ! includes all options for the nonlinear and linear solvers used
166: ! internally the timestepping routines.
167: !
169: call TSSetFromOptions(ts,ierr)
171: call TSSetUp(ts,ierr)
173: !
174: ! Perform the solve. This is where the timestepping takes place.
175: !
176: call TSSolve(ts,x,ierr)
177: call TSGetSolveTime(ts,ftime,ierr)
178: call TSGetStepNumber(ts,its,ierr)
179: call TSGetConvergedReason(ts,reason,ierr)
181: write(6,100) its,ftime,reason
182: 100 format('Number of pseudo time-steps ',i5,' final time ',1pe9.2 &
183: & ,' reason ',i3)
185: !
186: ! Free the data structures constructed above
187: !
189: call VecDestroy(x,ierr)
190: call VecDestroy(r,ierr)
191: call MatDestroy(J,ierr)
192: call TSDestroy(ts,ierr)
193: call PetscFinalize(ierr)
194: end
196: !
197: ! -------------------- Form initial approximation -----------------
198: !
199: subroutine FormInitialGuess(X,user,ierr)
200: use petscts
201: implicit none
203: Vec X
204: PetscReal user(3)
205: PetscInt i,j,row,mx,my
206: PetscErrorCode ierr
207: PetscOffset xidx
208: PetscReal one,lambda
209: PetscReal temp1,temp,hx,hy
210: PetscScalar xx(2)
211: PetscInt param,lmx,lmy
212: parameter (param = 1,lmx = 2,lmy = 3)
214: one = 1.0
216: mx = int(user(lmx))
217: my = int(user(lmy))
218: lambda = user(param)
220: hy = one / (my-1)
221: hx = one / (mx-1)
223: call VecGetArray(X,xx,xidx,ierr)
224: temp1 = lambda/(lambda + one)
225: do 10, j=1,my
226: temp = min(j-1,my-j)*hy
227: do 20 i=1,mx
228: row = i + (j-1)*mx
229: if (i .eq. 1 .or. j .eq. 1 .or. &
230: & i .eq. mx .or. j .eq. my) then
231: xx(row+xidx) = 0.0
232: else
233: xx(row+xidx) = &
234: & temp1*sqrt(min(min(i-1,mx-i)*hx,temp))
235: endif
236: 20 continue
237: 10 continue
238: call VecRestoreArray(X,xx,xidx,ierr)
239: return
240: end
241: !
242: ! -------------------- Evaluate Function F(x) ---------------------
243: !
244: subroutine FormFunction(ts,t,X,F,user,ierr)
245: use petscts
246: implicit none
248: TS ts
249: PetscReal t
250: Vec X,F
251: PetscReal user(3)
252: PetscErrorCode ierr
253: PetscInt i,j,row,mx,my
254: PetscOffset xidx,fidx
255: PetscReal two,lambda
256: PetscReal hx,hy,hxdhy,hydhx
257: PetscScalar ut,ub,ul,ur,u
258: PetscScalar uxx,uyy,sc
259: PetscScalar xx(2),ff(2)
260: PetscInt param,lmx,lmy
261: parameter (param = 1,lmx = 2,lmy = 3)
263: two = 2.0
265: mx = int(user(lmx))
266: my = int(user(lmy))
267: lambda = user(param)
269: hx = 1.0 / real(mx-1)
270: hy = 1.0 / real(my-1)
271: sc = hx*hy
272: hxdhy = hx/hy
273: hydhx = hy/hx
275: call VecGetArrayRead(X,xx,xidx,ierr)
276: call VecGetArray(F,ff,fidx,ierr)
277: do 10 j=1,my
278: do 20 i=1,mx
279: row = i + (j-1)*mx
280: if (i .eq. 1 .or. j .eq. 1 .or. &
281: & i .eq. mx .or. j .eq. my) then
282: ff(row+fidx) = xx(row+xidx)
283: else
284: u = xx(row + xidx)
285: ub = xx(row - mx + xidx)
286: ul = xx(row - 1 + xidx)
287: ut = xx(row + mx + xidx)
288: ur = xx(row + 1 + xidx)
289: uxx = (-ur + two*u - ul)*hydhx
290: uyy = (-ut + two*u - ub)*hxdhy
291: ff(row+fidx) = -uxx - uyy + sc*lambda*exp(u)
292: u = -uxx - uyy + sc*lambda*exp(u)
293: endif
294: 20 continue
295: 10 continue
297: call VecRestoreArrayRead(X,xx,xidx,ierr)
298: call VecRestoreArray(F,ff,fidx,ierr)
299: return
300: end
301: !
302: ! -------------------- Evaluate Jacobian of F(x) --------------------
303: !
304: subroutine FormJacobian(ts,ctime,X,JJ,B,user,ierr)
305: use petscts
306: implicit none
308: TS ts
309: Vec X
310: Mat JJ,B
311: PetscReal user(3),ctime
312: PetscErrorCode ierr
313: Mat jac
314: PetscOffset xidx
315: PetscInt i,j,row(1),mx,my
316: PetscInt col(5),i1,i5
317: PetscScalar two,one,lambda
318: PetscScalar v(5),sc,xx(2)
319: PetscReal hx,hy,hxdhy,hydhx
321: PetscInt param,lmx,lmy
322: parameter (param = 1,lmx = 2,lmy = 3)
324: i1 = 1
325: i5 = 5
326: jac = B
327: two = 2.0
328: one = 1.0
330: mx = int(user(lmx))
331: my = int(user(lmy))
332: lambda = user(param)
334: hx = 1.0 / real(mx-1)
335: hy = 1.0 / real(my-1)
336: sc = hx*hy
337: hxdhy = hx/hy
338: hydhx = hy/hx
340: call VecGetArrayRead(X,xx,xidx,ierr)
341: do 10 j=1,my
342: do 20 i=1,mx
343: !
344: ! When inserting into PETSc matrices, indices start at 0
345: !
346: row(1) = i - 1 + (j-1)*mx
347: if (i .eq. 1 .or. j .eq. 1 .or. &
348: & i .eq. mx .or. j .eq. my) then
349: call MatSetValues(jac,i1,row,i1,row,one,INSERT_VALUES,ierr)
350: else
351: v(1) = hxdhy
352: col(1) = row(1) - mx
353: v(2) = hydhx
354: col(2) = row(1) - 1
355: v(3) = -two*(hydhx+hxdhy)+sc*lambda*exp(xx(row(1)+1+xidx))
356: col(3) = row(1)
357: v(4) = hydhx
358: col(4) = row(1) + 1
359: v(5) = hxdhy
360: col(5) = row(1) + mx
361: call MatSetValues(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)
362: endif
363: 20 continue
364: 10 continue
365: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
366: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
367: call VecRestoreArray(X,xx,xidx,ierr)
368: return
369: end
371: !/*TEST
372: !
373: ! test:
374: ! TODO: broken
375: ! 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
376: !
377: !TEST*/