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