Actual source code: ex1f.F

petsc-master 2016-10-26
Report Typos and Errors
```  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:       if (ierr .ne. 0) then
69:          print*,'Unable to initialize PETSc'
70:          stop
71:       endif
72:       user(lmx)        = 4
73:       user(lmy)        = 4
74:       user(param)      = 6.0

76: !
77: !     Allow user to set the grid dimensions and nonlinearity parameter at run-time
78: !
79:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,                    &
80:      &                         '-mx',user(lmx),flg,ierr)
81:       call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,            &
82:      &                        '-my',user(lmy),flg,ierr)
83:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,           &
84:      &                        '-param',user(param),flg,ierr)
85:       if (user(param) .ge. param_max .or.                               &
86:      &                                user(param) .le. param_min) then
87:         print*,'Parameter is out of range'
88:       endif
89:       if (user(lmx) .gt. user(lmy)) then
90:         dt = .5/user(lmx)
91:       else
92:         dt = .5/user(lmy)
93:       endif
94:       call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,       &
95:      &                         '-dt',dt,flg,ierr)
96:       N          = int(user(lmx)*user(lmy))

98: !
99: !      Create vectors to hold the solution and function value
100: !
101:       call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr)
102:       call VecDuplicate(x,r,ierr)

104: !
105: !    Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
106: !    in the sparse matrix. Note that this is not the optimal strategy see
107: !    the Performance chapter of the users manual for information on
108: !    preallocating memory in sparse matrices.
109: !
110:       i5 = 5
111:       call MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,i5,PETSC_NULL_INTEGER,    &
112:      &     J,ierr)

114: !
115: !     Create timestepper context
116: !

118:       call TSCreate(PETSC_COMM_WORLD,ts,ierr)
119:       call TSSetProblemType(ts,TS_NONLINEAR,ierr)

121: !
122: !     Tell the timestepper context where to compute solutions
123: !

125:       call TSSetSolution(ts,x,ierr)

127: !
128: !    Provide the call-back for the nonlinear function we are
129: !    evaluating. Thus whenever the timestepping routines need the
130: !    function they will call this routine. Note the final argument
131: !    is the application context used by the call-back functions.
132: !

134:       call TSSetRHSFunction(ts,PETSC_NULL_OBJECT,FormFunction,user,ierr)

136: !
137: !     Set the Jacobian matrix and the function used to compute
138: !     Jacobians.
139: !

141:       call TSSetRHSJacobian(ts,J,J,FormJacobian,user,ierr)

143: !
144: !       For the initial guess for the problem
145: !

147:       call FormInitialGuess(x,user,ierr)

149: !
150: !       This indicates that we are using pseudo timestepping to
151: !     find a steady state solution to the nonlinear problem.
152: !

154:       call TSSetType(ts,TSPSEUDO,ierr)

156: !
157: !       Set the initial time to start at (this is arbitrary for
158: !     steady state problems and the initial timestep given above
159: !

161:       zero = 0.0
162:       call TSSetInitialTimeStep(ts,zero,dt,ierr)

164: !
165: !      Set a large number of timesteps and final duration time
166: !     to insure convergence to steady state.
167: !
168:       i1000 = 1000
169:       tmax  = 1.e12
170:       call TSSetDuration(ts,i1000,tmax,ierr)
171:       call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)

173: !
174: !      Set any additional options from the options database. This
175: !     includes all options for the nonlinear and linear solvers used
176: !     internally the timestepping routines.
177: !

179:       call TSSetFromOptions(ts,ierr)

181:       call TSSetUp(ts,ierr)

183: !
184: !      Perform the solve. This is where the timestepping takes place.
185: !
186:       call TSSolve(ts,x,ierr)
187:       call TSGetSolveTime(ts,ftime,ierr)
188:       call TSGetTimeStepNumber(ts,its,ierr)
189:       call TSGetConvergedReason(ts,reason,ierr)

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

195: !
196: !     Free the data structures constructed above
197: !

199:       call VecDestroy(x,ierr)
200:       call VecDestroy(r,ierr)
201:       call MatDestroy(J,ierr)
202:       call TSDestroy(ts,ierr)
203:       call PetscFinalize(ierr)
204:       end

206: !
207: !  --------------------  Form initial approximation -----------------
208: !
209:       subroutine FormInitialGuess(X,user,ierr)
210:       implicit none
211:  #include <petsc/finclude/petscsys.h>
212:  #include <petsc/finclude/petscvec.h>
213:  #include <petsc/finclude/petscmat.h>
214:  #include <petsc/finclude/petscpc.h>
215:  #include <petsc/finclude/petscts.h>
216:       Vec              X
217:       PetscReal user(3)
218:       PetscInt  i,j,row,mx,my
219:       PetscErrorCode ierr
220:       PetscOffset      xidx
221:       PetscReal one,lambda
222:       PetscReal temp1,temp,hx,hy
223:       PetscScalar      xx(2)
224:       PetscInt          param,lmx,lmy
225:       parameter (param = 1,lmx = 2,lmy = 3)

227:       one = 1.0

229:       mx     = int(user(lmx))
230:       my     = int(user(lmy))
231:       lambda = user(param)

233:       hy    = one / (my-1)
234:       hx    = one / (mx-1)

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

279:       two = 2.0

281:       mx     = int(user(lmx))
282:       my     = int(user(lmy))
283:       lambda = user(param)

285:       hx    = 1.0 / real(mx-1)
286:       hy    = 1.0 / real(my-1)
287:       sc    = hx*hy
288:       hxdhy = hx/hy
289:       hydhx = hy/hx

291:       call VecGetArrayRead(X,xx,xidx,ierr)
292:       call VecGetArray(F,ff,fidx,ierr)
293:       do 10 j=1,my
294:         do 20 i=1,mx
295:           row = i + (j-1)*mx
296:           if (i .eq. 1 .or. j .eq. 1 .or.                               &
297:      &        i .eq. mx .or. j .eq. my) then
298:             ff(row+fidx) = xx(row+xidx)
299:           else
300:             u            = xx(row + xidx)
301:             ub           = xx(row - mx + xidx)
302:             ul           = xx(row - 1 + xidx)
303:             ut           = xx(row + mx + xidx)
304:             ur           = xx(row + 1 + xidx)
305:             uxx          = (-ur + two*u - ul)*hydhx
306:             uyy          = (-ut + two*u - ub)*hxdhy
307:             ff(row+fidx) = -uxx - uyy + sc*lambda*exp(u)
308:             u =  -uxx - uyy + sc*lambda*exp(u)
309:          endif
310:  20   continue
311:  10   continue

313:       call VecRestoreArrayRead(X,xx,xidx,ierr)
314:       call VecRestoreArray(F,ff,fidx,ierr)
315:       return
316:       end
317: !
318: !  --------------------  Evaluate Jacobian of F(x) --------------------
319: !
320:       subroutine FormJacobian(ts,ctime,X,JJ,B,user,ierr)
321:       implicit none
322:  #include <petsc/finclude/petscsys.h>
323:  #include <petsc/finclude/petscvec.h>
324:  #include <petsc/finclude/petscmat.h>
325:  #include <petsc/finclude/petscpc.h>
326:  #include <petsc/finclude/petscts.h>
327:       TS               ts
328:       Vec              X
329:       Mat              JJ,B
330:       PetscReal user(3),ctime
331:       PetscErrorCode   ierr
332:       Mat              jac
333:       PetscOffset xidx
334:       PetscInt    i,j,row(1),mx,my
335:       PetscInt    col(5),i1,i5
336:       PetscScalar two,one,lambda
337:       PetscScalar v(5),sc,xx(2)
338:       PetscReal hx,hy,hxdhy,hydhx

340:       PetscInt  param,lmx,lmy
341:       parameter (param = 1,lmx = 2,lmy = 3)

343:       i1 = 1
344:       i5 = 5
345:       jac = B
346:       two = 2.0
347:       one = 1.0

349:       mx     = int(user(lmx))
350:       my     = int(user(lmy))
351:       lambda = user(param)

353:       hx    = 1.0 / real(mx-1)
354:       hy    = 1.0 / real(my-1)
355:       sc    = hx*hy
356:       hxdhy = hx/hy
357:       hydhx = hy/hx

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

```