Actual source code: ex1f.F
petsc-3.4.5 2014-06-29
1: !
2: ! Solves the time dependent Bratu problem using pseudo-timestepping
3: !
4: ! Concepts: TS^pseudo-timestepping
5: ! Concepts: pseudo-timestepping
6: ! Concepts: TS^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: ! "petscsys.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/petscsys.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: PetscBool flg
57: PetscErrorCode ierr
58: double precision param_max,param_min,dt,tmax,zero
59: double precision ftime
60: TSConvergedReason reason
62: i5 = 5
63: param_max = 6.81
64: param_min = 0
66: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
67: user(lmx) = 4
68: user(lmy) = 4
69: user(param) = 6.0
71: !
72: ! Allow user to set the grid dimensions and nonlinearity parameter at run-time
73: !
74: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-mx',user(lmx), &
75: & flg,ierr)
76: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-my',user(lmy), &
77: & flg,ierr)
78: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-param', &
79: & user(param),flg,ierr)
80: if (user(param) .ge. param_max .or. &
81: & user(param) .le. param_min) then
82: print*,'Parameter is out of range'
83: endif
84: if (user(lmx) .gt. user(lmy)) then
85: dt = .5/user(lmx)
86: else
87: dt = .5/user(lmy)
88: endif
89: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-dt',dt,flg,ierr)
90: N = int(user(lmx)*user(lmy))
92: !
93: ! Create vectors to hold the solution and function value
94: !
95: call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr)
96: call VecDuplicate(x,r,ierr)
98: !
99: ! Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
100: ! in the sparse matrix. Note that this is not the optimal strategy see
101: ! the Performance chapter of the users manual for information on
102: ! preallocating memory in sparse matrices.
103: !
104: i5 = 5
105: call MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,i5,PETSC_NULL_INTEGER, &
106: & J,ierr)
108: !
109: ! Create timestepper context
110: !
112: call TSCreate(PETSC_COMM_WORLD,ts,ierr)
113: call TSSetProblemType(ts,TS_NONLINEAR,ierr)
115: !
116: ! Tell the timestepper context where to compute solutions
117: !
119: call TSSetSolution(ts,x,ierr)
121: !
122: ! Provide the call-back for the nonlinear function we are
123: ! evaluating. Thus whenever the timestepping routines need the
124: ! function they will call this routine. Note the final argument
125: ! is the application context used by the call-back functions.
126: !
128: call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormFunction,user,ierr)
130: !
131: ! Set the Jacobian matrix and the function used to compute
132: ! Jacobians.
133: !
135: call TSSetRHSJacobian(ts,J,J,FormJacobian,user,ierr)
137: !
138: ! For the initial guess for the problem
139: !
141: call FormInitialGuess(x,user,ierr)
143: !
144: ! This indicates that we are using pseudo timestepping to
145: ! find a steady state solution to the nonlinear problem.
146: !
148: call TSSetType(ts,TSPSEUDO,ierr)
150: !
151: ! Set the initial time to start at (this is arbitrary for
152: ! steady state problems and the initial timestep given above
153: !
155: zero = 0.0
156: call TSSetInitialTimeStep(ts,zero,dt,ierr)
158: !
159: ! Set a large number of timesteps and final duration time
160: ! to insure convergence to steady state.
161: !
162: i1000 = 1000
163: tmax = 1.e12
164: call TSSetDuration(ts,i1000,tmax,ierr)
166: !
167: ! Set any additional options from the options database. This
168: ! includes all options for the nonlinear and linear solvers used
169: ! internally the the timestepping routines.
170: !
172: call TSSetFromOptions(ts,ierr)
174: call TSSetUp(ts,ierr)
176: !
177: ! Perform the solve. This is where the timestepping takes place.
178: !
179: call TSSolve(ts,x,ierr)
180: call TSGetSolveTime(ts,ftime,ierr)
181: call TSGetTimeStepNumber(ts,its,ierr)
182: call TSGetConvergedReason(ts,reason,ierr)
184: write(6,100) its,ftime,reason
185: 100 format('Number of pseudo time-steps ',i5,' final time ',1pe8.2 &
186: & ,' reason ',i3)
188: !
189: ! Free the data structures constructed above
190: !
192: call VecDestroy(x,ierr)
193: call VecDestroy(r,ierr)
194: call MatDestroy(J,ierr)
195: call TSDestroy(ts,ierr)
196: call PetscFinalize(ierr)
197: end
199: !
200: ! -------------------- Form initial approximation -----------------
201: !
202: subroutine FormInitialGuess(X,user,ierr)
203: implicit none
204: #include <finclude/petscsys.h>
205: #include <finclude/petscvec.h>
206: #include <finclude/petscmat.h>
207: #include <finclude/petscpc.h>
208: #include <finclude/petscts.h>
209: Vec X
210: double precision user(3)
211: PetscInt i,j,row,mx,my
212: PetscErrorCode ierr
213: PetscOffset xidx
214: double precision one,lambda
215: double precision temp1,temp,hx,hy
216: PetscScalar xx(1)
217: PetscInt param,lmx,lmy
218: parameter (param = 1,lmx = 2,lmy = 3)
220: one = 1.0
222: mx = int(user(lmx))
223: my = int(user(lmy))
224: lambda = user(param)
226: hy = one / (my-1)
227: hx = one / (mx-1)
229: call VecGetArray(X,xx,xidx,ierr)
230: temp1 = lambda/(lambda + one)
231: do 10, j=1,my
232: temp = dble(min(j-1,my-j))*hy
233: do 20 i=1,mx
234: row = i + (j-1)*mx
235: if (i .eq. 1 .or. j .eq. 1 .or. &
236: & i .eq. mx .or. j .eq. my) then
237: xx(row+xidx) = 0.0
238: else
239: xx(row+xidx) = &
240: & temp1*sqrt(min(dble(min(i-1,mx-i))*hx,temp))
241: endif
242: 20 continue
243: 10 continue
244: call VecRestoreArray(X,xx,xidx,ierr)
245: return
246: end
247: !
248: ! -------------------- Evaluate Function F(x) ---------------------
249: !
250: subroutine FormFunction(ts,t,X,F,user,ierr)
251: implicit none
252: #include <finclude/petscsys.h>
253: #include <finclude/petscvec.h>
254: #include <finclude/petscmat.h>
255: #include <finclude/petscpc.h>
256: #include <finclude/petscts.h>
257: TS ts
258: double precision t
259: Vec X,F
260: double precision user(3)
261: PetscErrorCode ierr
262: PetscInt i,j,row,mx,my
263: PetscOffset xidx,fidx
264: double precision two,lambda
265: double precision hx,hy,hxdhy,hydhx
266: PetscScalar ut,ub,ul,ur,u
267: PetscScalar uxx,uyy,sc
268: PetscScalar xx(1),ff(1)
269: PetscInt param,lmx,lmy
270: parameter (param = 1,lmx = 2,lmy = 3)
272: two = 2.0
274: mx = int(user(lmx))
275: my = int(user(lmy))
276: lambda = user(param)
278: hx = 1.0 / dble(mx-1)
279: hy = 1.0 / dble(my-1)
280: sc = hx*hy
281: hxdhy = hx/hy
282: hydhx = hy/hx
284: call VecGetArray(X,xx,xidx,ierr)
285: call VecGetArray(F,ff,fidx,ierr)
286: do 10 j=1,my
287: do 20 i=1,mx
288: row = i + (j-1)*mx
289: if (i .eq. 1 .or. j .eq. 1 .or. &
290: & i .eq. mx .or. j .eq. my) then
291: ff(row+fidx) = xx(row+xidx)
292: else
293: u = xx(row + xidx)
294: ub = xx(row - mx + xidx)
295: ul = xx(row - 1 + xidx)
296: ut = xx(row + mx + xidx)
297: ur = xx(row + 1 + xidx)
298: uxx = (-ur + two*u - ul)*hydhx
299: uyy = (-ut + two*u - ub)*hxdhy
300: ff(row+fidx) = -uxx - uyy + sc*lambda*exp(u)
301: u = -uxx - uyy + sc*lambda*exp(u)
302: endif
303: 20 continue
304: 10 continue
306: call VecRestoreArray(X,xx,xidx,ierr)
307: call VecRestoreArray(F,ff,fidx,ierr)
308: return
309: end
310: !
311: ! -------------------- Evaluate Jacobian of F(x) --------------------
312: !
313: subroutine FormJacobian(ts,ctime,X,JJ,B,flag,user,ierr)
314: implicit none
315: #include <finclude/petscsys.h>
316: #include <finclude/petscvec.h>
317: #include <finclude/petscmat.h>
318: #include <finclude/petscpc.h>
319: #include <finclude/petscts.h>
320: TS ts
321: Vec X
322: Mat JJ,B
323: MatStructure flag
324: double precision user(3),ctime
325: PetscErrorCode ierr
326: Mat jac
327: PetscOffset xidx
328: PetscInt i,j,row(1),mx,my
329: PetscInt col(5),i1,i5
330: PetscScalar two,one,lambda
331: PetscScalar v(5),sc,xx(1)
332: double precision hx,hy,hxdhy,hydhx
334: PetscInt param,lmx,lmy
335: parameter (param = 1,lmx = 2,lmy = 3)
337: i1 = 1
338: i5 = 5
339: jac = B
340: two = 2.0
341: one = 1.0
343: mx = int(user(lmx))
344: my = int(user(lmy))
345: lambda = user(param)
347: hx = 1.0 / dble(mx-1)
348: hy = 1.0 / dble(my-1)
349: sc = hx*hy
350: hxdhy = hx/hy
351: hydhx = hy/hx
353: call VecGetArray(X,xx,xidx,ierr)
354: do 10 j=1,my
355: do 20 i=1,mx
356: !
357: ! When inserting into PETSc matrices, indices start at 0
358: !
359: row(1) = i - 1 + (j-1)*mx
360: if (i .eq. 1 .or. j .eq. 1 .or. &
361: & i .eq. mx .or. j .eq. my) then
362: call MatSetValues(jac,i1,row,i1,row,one,INSERT_VALUES,ierr)
363: else
364: v(1) = hxdhy
365: col(1) = row(1) - mx
366: v(2) = hydhx
367: col(2) = row(1) - 1
368: v(3) = -two*(hydhx+hxdhy)+sc*lambda*exp(xx(row(1)+1+xidx))
369: col(3) = row(1)
370: v(4) = hydhx
371: col(4) = row(1) + 1
372: v(5) = hxdhy
373: col(5) = row(1) + mx
374: call MatSetValues(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)
375: endif
376: 20 continue
377: 10 continue
378: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
379: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
380: call VecRestoreArray(X,xx,xidx,ierr)
381: flag = SAME_NONZERO_PATTERN
382: return
383: end