Actual source code: ex1f.F

petsc-3.3-p7 2013-05-11
```  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: !           "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
70:
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))
91:
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,ftime,ierr)
180:       call TSGetTimeStepNumber(ts,its,ierr)
181:       call TSGetConvergedReason(ts,reason,ierr)

183:       write(6,100) its,ftime,reason
184:  100  format('Number of pseudo time-steps ',i5,' final time ',1pe8.2    &
185:      &     ,' reason ',i3)

187: !
188: !     Free the data structures constructed above
189: !

191:       call VecDestroy(x,ierr)
192:       call VecDestroy(r,ierr)
193:       call MatDestroy(J,ierr)
194:       call TSDestroy(ts,ierr)
195:       call PetscFinalize(ierr)
196:       end

198: !
199: !  --------------------  Form initial approximation -----------------
200: !
201:       subroutine FormInitialGuess(X,user,ierr)
202:       implicit none
203: #include <finclude/petscsys.h>
204: #include <finclude/petscvec.h>
205: #include <finclude/petscmat.h>
206: #include <finclude/petscpc.h>
207: #include <finclude/petscts.h>
208:       Vec              X
209:       double precision user(3)
210:       PetscInt  i,j,row,mx,my
211:       PetscErrorCode ierr
212:       PetscOffset      xidx
213:       double precision one,lambda
214:       double precision temp1,temp,hx,hy
215:       PetscScalar      xx(1)
216:       PetscInt          param,lmx,lmy
217:       parameter (param = 1,lmx = 2,lmy = 3)

219:       one = 1.0

221:       mx     = int(user(lmx))
222:       my     = int(user(lmy))
223:       lambda = user(param)

225:       hy    = one / (my-1)
226:       hx    = one / (mx-1)

228:       call VecGetArray(X,xx,xidx,ierr)
229:       temp1 = lambda/(lambda + one)
230:       do 10, j=1,my
231:         temp = dble(min(j-1,my-j))*hy
232:         do 20 i=1,mx
233:           row = i + (j-1)*mx
234:           if (i .eq. 1 .or. j .eq. 1 .or.                               &
235:      &        i .eq. mx .or. j .eq. my) then
236:             xx(row+xidx) = 0.0
237:           else
238:             xx(row+xidx) =                                              &
239:      &        temp1*sqrt(min(dble(min(i-1,mx-i))*hx,temp))
240:           endif
241:  20     continue
242:  10   continue
243:       call VecRestoreArray(X,xx,xidx,ierr)
244:       return
245:       end
246: !
247: !  --------------------  Evaluate Function F(x) ---------------------
248: !
249:       subroutine FormFunction(ts,t,X,F,user,ierr)
250:       implicit none
251: #include <finclude/petscsys.h>
252: #include <finclude/petscvec.h>
253: #include <finclude/petscmat.h>
254: #include <finclude/petscpc.h>
255: #include <finclude/petscts.h>
256:       TS       ts
257:       double precision  t
258:       Vec               X,F
259:       double precision  user(3)
260:       PetscErrorCode     ierr
261:       PetscInt         i,j,row,mx,my
262:       PetscOffset       xidx,fidx
263:       double precision  two,lambda
264:       double precision  hx,hy,hxdhy,hydhx
265:       PetscScalar  ut,ub,ul,ur,u
266:       PetscScalar  uxx,uyy,sc
267:       PetscScalar  xx(1),ff(1)
268:       PetscInt     param,lmx,lmy
269:       parameter (param = 1,lmx = 2,lmy = 3)

271:       two = 2.0

273:       mx     = int(user(lmx))
274:       my     = int(user(lmy))
275:       lambda = user(param)

277:       hx    = 1.0 / dble(mx-1)
278:       hy    = 1.0 / dble(my-1)
279:       sc    = hx*hy
280:       hxdhy = hx/hy
281:       hydhx = hy/hx

283:       call VecGetArray(X,xx,xidx,ierr)
284:       call VecGetArray(F,ff,fidx,ierr)
285:       do 10 j=1,my
286:         do 20 i=1,mx
287:           row = i + (j-1)*mx
288:           if (i .eq. 1 .or. j .eq. 1 .or.                               &
289:      &        i .eq. mx .or. j .eq. my) then
290:             ff(row+fidx) = xx(row+xidx)
291:           else
292:             u            = xx(row + xidx)
293:             ub           = xx(row - mx + xidx)
294:             ul           = xx(row - 1 + xidx)
295:             ut           = xx(row + mx + xidx)
296:             ur           = xx(row + 1 + xidx)
297:             uxx          = (-ur + two*u - ul)*hydhx
298:             uyy          = (-ut + two*u - ub)*hxdhy
299:             ff(row+fidx) = -uxx - uyy + sc*lambda*exp(u)
300:             u =  -uxx - uyy + sc*lambda*exp(u)
301:          endif
302:  20   continue
303:  10   continue

305:       call VecRestoreArray(X,xx,xidx,ierr)
306:       call VecRestoreArray(F,ff,fidx,ierr)
307:       return
308:       end
309: !
310: !  --------------------  Evaluate Jacobian of F(x) --------------------
311: !
312:       subroutine FormJacobian(ts,ctime,X,JJ,B,flag,user,ierr)
313:       implicit none
314: #include <finclude/petscsys.h>
315: #include <finclude/petscvec.h>
316: #include <finclude/petscmat.h>
317: #include <finclude/petscpc.h>
318: #include <finclude/petscts.h>
319:       TS               ts
320:       Vec              X
321:       Mat              JJ,B
322:       MatStructure     flag
323:       double precision user(3),ctime
324:       PetscErrorCode   ierr
325:       Mat              jac
326:       PetscOffset xidx
327:       PetscInt    i,j,row(1),mx,my
328:       PetscInt    col(5),i1,i5
329:       PetscScalar two,one,lambda
330:       PetscScalar v(5),sc,xx(1)
331:       double precision hx,hy,hxdhy,hydhx

333:       PetscInt  param,lmx,lmy
334:       parameter (param = 1,lmx = 2,lmy = 3)

336:       i1 = 1
337:       i5 = 5
338:       jac = B
339:       two = 2.0
340:       one = 1.0

342:       mx     = int(user(lmx))
343:       my     = int(user(lmy))
344:       lambda = user(param)

346:       hx    = 1.0 / dble(mx-1)
347:       hy    = 1.0 / dble(my-1)
348:       sc    = hx*hy
349:       hxdhy = hx/hy
350:       hydhx = hy/hx

352:       call VecGetArray(X,xx,xidx,ierr)
353:       do 10 j=1,my
354:         do 20 i=1,mx
355: !
356: !      When inserting into PETSc matrices, indices start at 0
357: !
358:           row(1) = i - 1 + (j-1)*mx
359:           if (i .eq. 1 .or. j .eq. 1 .or.                               &
360:      &        i .eq. mx .or. j .eq. my) then
361:             call MatSetValues(jac,i1,row,i1,row,one,INSERT_VALUES,ierr)
362:           else
363:             v(1)   = hxdhy
364:             col(1) = row(1) - mx
365:             v(2)   = hydhx
366:             col(2) = row(1) - 1
367:             v(3)   = -two*(hydhx+hxdhy)+sc*lambda*exp(xx(row(1)+1+xidx))
368:             col(3) = row(1)
369:             v(4)   = hydhx
370:             col(4) = row(1) + 1
371:             v(5)   = hxdhy
372:             col(5) = row(1) + mx
373:             call MatSetValues(jac,i1,row,i5,col,v,INSERT_VALUES,ierr)
374:           endif
375:  20     continue
376:  10   continue
377:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
378:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
379:       call VecRestoreArray(X,xx,xidx,ierr)
380:       flag = SAME_NONZERO_PATTERN
381:       return
382:       end

```