Actual source code: ex1f.F

petsc-master 2017-04-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: !
 23: !23456789012345678901234567890123456789012345678901234567890123456789012
 24:       program main
 25:  #include <petsc/finclude/petscts.h>
 26:       use petscts
 27:       implicit none

 29: !
 30: !  Create an application context to contain data needed by the
 31: !  application-provided call-back routines, FormJacobian() and
 32: !  FormFunction(). We use a double precision array with three
 33: !  entries indexed by param, lmx, lmy.
 34: !
 35:       PetscReal user(3)
 36:       PetscInt          param,lmx,lmy,i5
 37:       parameter (param = 1,lmx = 2,lmy = 3)
 38: !
 39: !   User-defined routines
 40: !
 41:       external FormJacobian,FormFunction
 42: !
 43: !   Data for problem
 44: !
 45:       TS                ts
 46:       Vec               x,r
 47:       Mat               J
 48:       PetscInt           its,N,i1000,itmp
 49:       PetscBool  flg
 50:       PetscErrorCode      ierr
 51:       PetscReal  param_max,param_min,dt
 52:       PetscReal  tmax,zero
 53:       PetscReal  ftime
 54:       TSConvergedReason reason

 56:       i5 = 5
 57:       param_max = 6.81
 58:       param_min = 0

 60:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 61:       if (ierr .ne. 0) then
 62:          print*,'Unable to initialize PETSc'
 63:          stop
 64:       endif
 65:       user(lmx)        = 4
 66:       user(lmy)        = 4
 67:       user(param)      = 6.0

 69: !
 70: !     Allow user to set the grid dimensions and nonlinearity parameter at run-time
 71: !
 72:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,                    &
 73:      &     '-mx',user(lmx),flg,ierr)
 74:       itmp = 4
 75:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,            &
 76:      &                        '-my',itmp,flg,ierr)
 77:       user(lmy) = itmp
 78:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,           &
 79:      &                        '-param',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_OPTIONS,PETSC_NULL_CHARACTER,       &
 90:      &                         '-dt',dt,flg,ierr)
 91:       N          = int(user(lmx)*user(lmy))

 93: !
 94: !      Create vectors to hold the solution and function value
 95: !
 96:       call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr)
 97:       call VecDuplicate(x,r,ierr)

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

109: !
110: !     Create timestepper context
111: !

113:       call TSCreate(PETSC_COMM_WORLD,ts,ierr)
114:       call TSSetProblemType(ts,TS_NONLINEAR,ierr)

116: !
117: !     Tell the timestepper context where to compute solutions
118: !

120:       call TSSetSolution(ts,x,ierr)

122: !
123: !    Provide the call-back for the nonlinear function we are
124: !    evaluating. Thus whenever the timestepping routines need the
125: !    function they will call this routine. Note the final argument
126: !    is the application context used by the call-back functions.
127: !

129:       call TSSetRHSFunction(ts,PETSC_NULL_VEC,FormFunction,user,ierr)

131: !
132: !     Set the Jacobian matrix and the function used to compute
133: !     Jacobians.
134: !

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

138: !
139: !       For the initial guess for the problem
140: !

142:       call FormInitialGuess(x,user,ierr)

144: !
145: !       This indicates that we are using pseudo timestepping to
146: !     find a steady state solution to the nonlinear problem.
147: !

149:       call TSSetType(ts,TSPSEUDO,ierr)

151: !
152: !       Set the initial time to start at (this is arbitrary for
153: !     steady state problems and the initial timestep given above
154: !

156:       zero = 0.0
157:       call TSSetInitialTimeStep(ts,zero,dt,ierr)

159: !
160: !      Set a large number of timesteps and final duration time
161: !     to insure convergence to steady state.
162: !
163:       i1000 = 1000
164:       tmax  = 1.e12
165:       call TSSetDuration(ts,i1000,tmax,ierr)
166:       call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)

168: !
169: !      Set any additional options from the options database. This
170: !     includes all options for the nonlinear and linear solvers used
171: !     internally the timestepping routines.
172: !

174:       call TSSetFromOptions(ts,ierr)

176:       call TSSetUp(ts,ierr)

178: !
179: !      Perform the solve. This is where the timestepping takes place.
180: !
181:       call TSSolve(ts,x,ierr)
182:       call TSGetSolveTime(ts,ftime,ierr)
183:       call TSGetTimeStepNumber(ts,its,ierr)
184:       call TSGetConvergedReason(ts,reason,ierr)

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

190: !
191: !     Free the data structures constructed above
192: !

194:       call VecDestroy(x,ierr)
195:       call VecDestroy(r,ierr)
196:       call MatDestroy(J,ierr)
197:       call TSDestroy(ts,ierr)
198:       call PetscFinalize(ierr)
199:       end

201: !
202: !  --------------------  Form initial approximation -----------------
203: !
204:       subroutine FormInitialGuess(X,user,ierr)
205:       use petscts
206:       implicit none

208:       Vec              X
209:       PetscReal user(3)
210:       PetscInt  i,j,row,mx,my
211:       PetscErrorCode ierr
212:       PetscOffset      xidx
213:       PetscReal one,lambda
214:       PetscReal temp1,temp,hx,hy
215:       PetscScalar      xx(2)
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 = 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(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:       use petscts
251:       implicit none

253:       TS       ts
254:       PetscReal  t
255:       Vec               X,F
256:       PetscReal  user(3)
257:       PetscErrorCode     ierr
258:       PetscInt         i,j,row,mx,my
259:       PetscOffset       xidx,fidx
260:       PetscReal  two,lambda
261:       PetscReal  hx,hy,hxdhy,hydhx
262:       PetscScalar  ut,ub,ul,ur,u
263:       PetscScalar  uxx,uyy,sc
264:       PetscScalar  xx(2),ff(2)
265:       PetscInt     param,lmx,lmy
266:       parameter (param = 1,lmx = 2,lmy = 3)

268:       two = 2.0

270:       mx     = int(user(lmx))
271:       my     = int(user(lmy))
272:       lambda = user(param)

274:       hx    = 1.0 / real(mx-1)
275:       hy    = 1.0 / real(my-1)
276:       sc    = hx*hy
277:       hxdhy = hx/hy
278:       hydhx = hy/hx

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

302:       call VecRestoreArrayRead(X,xx,xidx,ierr)
303:       call VecRestoreArray(F,ff,fidx,ierr)
304:       return
305:       end
306: !
307: !  --------------------  Evaluate Jacobian of F(x) --------------------
308: !
309:       subroutine FormJacobian(ts,ctime,X,JJ,B,user,ierr)
310:       use petscts
311:       implicit none

313:       TS               ts
314:       Vec              X
315:       Mat              JJ,B
316:       PetscReal user(3),ctime
317:       PetscErrorCode   ierr
318:       Mat              jac
319:       PetscOffset xidx
320:       PetscInt    i,j,row(1),mx,my
321:       PetscInt    col(5),i1,i5
322:       PetscScalar two,one,lambda
323:       PetscScalar v(5),sc,xx(2)
324:       PetscReal hx,hy,hxdhy,hydhx

326:       PetscInt  param,lmx,lmy
327:       parameter (param = 1,lmx = 2,lmy = 3)

329:       i1 = 1
330:       i5 = 5
331:       jac = B
332:       two = 2.0
333:       one = 1.0

335:       mx     = int(user(lmx))
336:       my     = int(user(lmy))
337:       lambda = user(param)

339:       hx    = 1.0 / real(mx-1)
340:       hy    = 1.0 / real(my-1)
341:       sc    = hx*hy
342:       hxdhy = hx/hy
343:       hydhx = hy/hx

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