Actual source code: ex1f.F

petsc-master 2016-08-29
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:       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(1)
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(1),ff(1)
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(1)
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