Actual source code: ex1f.F
1: !
2: ! Solves the time dependent Bratu problem using pseudo-timestepping
3: !
4: ! Concepts: TS^pseudo-timestepping
5: ! Concepts: pseudo-timestepping
6: ! Concepts: nonlinear problems
7: ! Processors: 1
8: !
9: ! This code demonstrates how one may solve a nonlinear problem
10: ! with pseudo-timestepping. In this simple example, the pseudo-timestep
11: ! is the same for all grid points, i.e., this is equivalent to using
12: ! the backward Euler method with a variable timestep.
13: !
14: ! Note: This example does not require pseudo-timestepping since it
15: ! is an easy nonlinear problem, but it is included to demonstrate how
16: ! the pseudo-timestepping may be done.
17: !
18: ! See snes/examples/tutorials/ex4.c[ex4f.F] and
19: ! snes/examples/tutorials/ex5.c[ex5f.F] where the problem is described
20: ! and solved using the method of Newton alone.
21: !
22: ! Include "petscts.h" to use the PETSc timestepping routines,
23: ! "petsc.h" for basic PETSc operation,
24: ! "petscmat.h" for matrix operations,
25: ! "petscpc.h" for preconditions, and
26: ! "petscvec.h" for vector operations.
27: !
28: !23456789012345678901234567890123456789012345678901234567890123456789012
29: program main
30: implicit none
31: #include finclude/petsc.h
32: #include finclude/petscvec.h
33: #include finclude/petscmat.h
34: #include finclude/petscpc.h
35: #include finclude/petscts.h
36: !
37: ! Create an application context to contain data needed by the
38: ! application-provided call-back routines, FormJacobian() and
39: ! FormFunction(). We use a double precision array with three
40: ! entries indexed by param, lmx, lmy.
41: !
42: double precision user(3)
43: PetscInt param,lmx,lmy,i5
44: parameter (param = 1,lmx = 2,lmy = 3)
45: !
46: ! User-defined routines
47: !
48: external FormJacobian,FormFunction
49: !
50: ! Data for problem
51: !
52: TS ts
53: Vec x,r
54: Mat J
55: PetscInt its,N,i1000
56: PetscTruth flg
57: PetscErrorCode ierr
58: double precision param_max,param_min,dt,tmax,zero
59: double precision ftime
61: i5 = 5
62: param_max = 6.81
63: param_min = 0
65: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
66: user(lmx) = 4
67: user(lmy) = 4
68: user(param) = 6.0
69:
70: !
71: ! Allow user to set the grid dimensions and nonlinearity parameter at run-time
72: !
73: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-mx',user(lmx), &
74: & flg,ierr)
75: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-my',user(lmy), &
76: & flg,ierr)
77: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-param', &
78: & user(param),flg,ierr)
79: if (user(param) .ge. param_max .or. &
80: & user(param) .le. param_min) then
81: print*,'Parameter is out of range'
82: endif
83: if (user(lmx) .gt. user(lmy)) then
84: dt = .5/user(lmx)
85: else
86: dt = .5/user(lmy)
87: endif
88: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)
89: N = int(user(lmx)*user(lmy))
90:
91: !
92: ! Create vectors to hold the solution and function value
93: !
94: call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr)
95: call VecDuplicate(x,r,ierr)
97: !
98: ! Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
99: ! in the sparse matrix. Note that this is not the optimal strategy see
100: ! the Performance chapter of the users manual for information on
101: ! preallocating memory in sparse matrices.
102: !
103: i5 = 5
104: call MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,i5,PETSC_NULL_INTEGER, &
105: & J,ierr)
107: !
108: ! Create timestepper context
109: !
111: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
112: call TSSetProblemType(ts,TS_NONLINEAR,ierr)
114: !
115: ! Tell the timestepper context where to compute solutions
116: !
118: call TSSetSolution(ts,x,ierr)
120: !
121: ! Provide the call-back for the nonlinear function we are
122: ! evaluating. Thus whenever the timestepping routines need the
123: ! function they will call this routine. Note the final argument
124: ! is the application context used by the call-back functions.
125: !
127: call TSSetRHSFunction(ts,FormFunction,user,ierr)
129: !
130: ! Set the Jacobian matrix and the function used to compute
131: ! Jacobians.
132: !
134: call TSSetRHSJacobian(ts,J,J,FormJacobian,user,ierr)
136: !
137: ! For the initial guess for the problem
138: !
140: call FormInitialGuess(x,user,ierr)
142: !
143: ! This indicates that we are using pseudo timestepping to
144: ! find a steady state solution to the nonlinear problem.
145: !
147: call TSSetType(ts,TS_PSEUDO,ierr)
149: !
150: ! Set the initial time to start at (this is arbitrary for
151: ! steady state problems and the initial timestep given above
152: !
154: zero = 0.0
155: call TSSetInitialTimeStep(ts,zero,dt,ierr)
157: !
158: ! Set a large number of timesteps and final duration time
159: ! to insure convergence to steady state.
160: !
161: i1000 = 1000
162: tmax = 1.e12
163: call TSSetDuration(ts,i1000,tmax,ierr)
165: !
166: ! Set any additional options from the options database. This
167: ! includes all options for the nonlinear and linear solvers used
168: ! internally the the timestepping routines.
169: !
171: call TSSetFromOptions(ts,ierr)
173: call TSSetUp(ts,ierr)
175: !
176: ! Perform the solve. This is where the timestepping takes place.
177: !
178:
179: call TSStep(ts,its,ftime,ierr)
180:
181: write(6,100) its,ftime
182: 100 format('Number of pseudo time-steps ',i5,' final time ',1pe8.2)
184: !
185: ! Free the data structures constructed above
186: !
188: call VecDestroy(x,ierr)
189: call VecDestroy(r,ierr)
190: call MatDestroy(J,ierr)
191: call TSDestroy(ts,ierr)
192: call PetscFinalize(ierr)
193: end
195: !
196: ! -------------------- Form initial approximation -----------------
197: !
198: subroutine FormInitialGuess(X,user,ierr)
199: implicit none
200: #include finclude/petsc.h
201: #include finclude/petscvec.h
202: #include finclude/petscmat.h
203: #include finclude/petscpc.h
204: #include finclude/petscts.h
205: Vec X
206: double precision user(3)
207: PetscInt i,j,row,mx,my
208: PetscErrorCode ierr
209: PetscOffset xidx
210: double precision one,lambda
211: double precision temp1,temp,hx,hy
212: PetscScalar xx(1)
213: PetscInt param,lmx,lmy
214: parameter (param = 1,lmx = 2,lmy = 3)
216: one = 1.0
218: mx = int(user(lmx))
219: my = int(user(lmy))
220: lambda = user(param)
222: hy = one / (my-1)
223: hx = one / (mx-1)
225: call VecGetArray(X,xx,xidx,ierr)
226: temp1 = lambda/(lambda + one)
227: do 10, j=1,my
228: temp = dble(min(j-1,my-j))*hy
229: do 20 i=1,mx
230: row = i + (j-1)*mx
231: if (i .eq. 1 .or. j .eq. 1 .or. &
232: & i .eq. mx .or. j .eq. my) then
233: xx(row+xidx) = 0.0
234: else
235: xx(row+xidx) = &
236: & temp1*sqrt(min(dble(min(i-1,mx-i))*hx,temp))
237: endif
238: 20 continue
239: 10 continue
240: call VecRestoreArray(X,xx,xidx,ierr)
241: return
242: end
243: !
244: ! -------------------- Evaluate Function F(x) ---------------------
245: !
246: subroutine FormFunction(ts,t,X,F,user,ierr)
247: implicit none
248: #include finclude/petsc.h
249: #include finclude/petscvec.h
250: #include finclude/petscmat.h
251: #include finclude/petscpc.h
252: #include finclude/petscts.h
253: TS ts
254: double precision t
255: Vec X,F
256: double precision user(3)
257: PetscErrorCode ierr
258: PetscInt i,j,row,mx,my
259: PetscOffset xidx,fidx
260: double precision two,lambda
261: double precision hx,hy,hxdhy,hydhx
262: PetscScalar ut,ub,ul,ur,u
263: PetscScalar uxx,uyy,sc
264: PetscScalar xx(1),ff(1)
265: PetscInt param,lmx,lmy
266: parameter (param = 1,lmx = 2,lmy = 3)
268: two = 2.0
270: mx = int(user(lmx))
271: my = int(user(lmy))
272: lambda = user(param)
274: hx = 1.0 / dble(mx-1)
275: hy = 1.0 / dble(my-1)
276: sc = hx*hy
277: hxdhy = hx/hy
278: hydhx = hy/hx
280: call VecGetArray(X,xx,xidx,ierr)
281: call VecGetArray(F,ff,fidx,ierr)
282: do 10 j=1,my
283: do 20 i=1,mx
284: row = i + (j-1)*mx
285: if (i .eq. 1 .or. j .eq. 1 .or. &
286: & i .eq. mx .or. j .eq. my) then
287: ff(row+fidx) = xx(row+xidx)
288: else
289: u = xx(row + xidx)
290: ub = xx(row - mx + xidx)
291: ul = xx(row - 1 + xidx)
292: ut = xx(row + mx + xidx)
293: ur = xx(row + 1 + xidx)
294: uxx = (-ur + two*u - ul)*hydhx
295: uyy = (-ut + two*u - ub)*hxdhy
296: ff(row+fidx) = -uxx - uyy + sc*lambda*exp(u)
297: u = -uxx - uyy + sc*lambda*exp(u)
298: endif
299: 20 continue
300: 10 continue
302: call VecRestoreArray(X,xx,xidx,ierr)
303: call VecRestoreArray(F,ff,fidx,ierr)
304: return
305: end
306: !
307: ! -------------------- Evaluate Jacobian of F(x) --------------------
308: !
309: subroutine FormJacobian(ts,ctime,X,JJ,B,flag,user,ierr)
310: implicit none
311: #include finclude/petsc.h
312: #include finclude/petscvec.h
313: #include finclude/petscmat.h
314: #include finclude/petscpc.h
315: #include finclude/petscts.h
316: TS ts
317: Vec X
318: Mat JJ,B
319: MatStructure flag
320: double precision user(3),ctime
321: PetscErrorCode ierr
322: Mat jac
323: PetscOffset xidx
324: PetscInt i,j,row(1),mx,my
325: PetscInt col(5),i1,i5
326: PetscScalar two,one,lambda
327: PetscScalar v(5),sc,xx(1)
328: double precision hx,hy,hxdhy,hydhx
330: PetscInt param,lmx,lmy
331: parameter (param = 1,lmx = 2,lmy = 3)
333: i1 = 1
334: i5 = 5
335: jac = B
336: two = 2.0
337: one = 1.0
339: mx = int(user(lmx))
340: my = int(user(lmy))
341: lambda = user(param)
343: hx = 1.0 / dble(mx-1)
344: hy = 1.0 / dble(my-1)
345: sc = hx*hy
346: hxdhy = hx/hy
347: hydhx = hy/hx
349: call VecGetArray(X,xx,xidx,ierr)
350: do 10 j=1,my
351: do 20 i=1,mx
352: !
353: ! When inserting into PETSc matrices, indices start at 0
354: !
355: row(1) = i - 1 + (j-1)*mx
356: if (i .eq. 1 .or. j .eq. 1 .or. &
357: & i .eq. mx .or. j .eq. my) then
358: call MatSetValues(jac,i1,row,i1,row,one,INSERT_VALUES,ierr)
359: else
360: v(1) = hxdhy
361: col(1) = row(1) - mx
362: v(2) = hydhx
363: col(2) = row(1) - 1
364: v(3) = -two*(hydhx+hxdhy)+sc*lambda*exp(xx(row(1)+1+xidx))
365: col(3) = row(1)
366: v(4) = hydhx
367: col(4) = row(1) + 1
368: v(5) = hxdhy
369: col(5) = row(1) + mx
370: call MatSetValues(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)
371: endif
372: 20 continue
373: 10 continue
374: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
375: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
376: call VecRestoreArray(X,xx,xidx,ierr)
377: flag = SAME_NONZERO_PATTERN
378: return
379: end