Actual source code: ex1f.F

  1: !
  2: !   Solves the time dependent Bratu problem using pseudo-timestepping
  3: !
  4: !   This code demonstrates how one may solve a nonlinear problem
  5: !   with pseudo-timestepping. In this simple example, the pseudo-timestep
  6: !   is the same for all grid points, i.e., this is equivalent to using
  7: !   the backward Euler method with a variable timestep.
  8: !
  9: !   Note: This example does not require pseudo-timestepping since it
 10: !   is an easy nonlinear problem, but it is included to demonstrate how
 11: !   the pseudo-timestepping may be done.
 12: !
 13: !   See snes/tutorials/ex4.c[ex4f.F] and
 14: !   snes/tutorials/ex5.c[ex5f.F] where the problem is described
 15: !   and solved using the method of Newton alone.
 16: !
 17: !
 18: !23456789012345678901234567890123456789012345678901234567890123456789012
 19:       program main
 20: #include <petsc/finclude/petscts.h>
 21:       use petscts
 22:       implicit none

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

 51:       i5 = 5
 52:       param_max = 6.81
 53:       param_min = 0

 55:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 56:       if (ierr .ne. 0) then
 57:          print*,'Unable to initialize PETSc'
 58:          stop
 59:       endif
 60:       user(lmx)        = 4
 61:       user(lmy)        = 4
 62:       user(param)      = 6.0

 64: !
 65: !     Allow user to set the grid dimensions and nonlinearity parameter at run-time
 66: !
 67:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,                    &
 68:      &     '-mx',user(lmx),flg,ierr)
 69:       itmp = 4
 70:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,            &
 71:      &                        '-my',itmp,flg,ierr)
 72:       user(lmy) = itmp
 73:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,           &
 74:      &                        '-param',user(param),flg,ierr)
 75:       if (user(param) .ge. param_max .or.                               &
 76:      &                                user(param) .le. param_min) then
 77:         print*,'Parameter is out of range'
 78:       endif
 79:       if (user(lmx) .gt. user(lmy)) then
 80:         dt = .5/user(lmx)
 81:       else
 82:         dt = .5/user(lmy)
 83:       endif
 84:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,       &
 85:      &                         '-dt',dt,flg,ierr)
 86:       N          = int(user(lmx)*user(lmy))

 88: !
 89: !      Create vectors to hold the solution and function value
 90: !
 91:       call VecCreateSeq(PETSC_COMM_SELF,N,x,ierr)
 92:       call VecDuplicate(x,r,ierr)

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

104: !
105: !     Create timestepper context
106: !

108:       call TSCreate(PETSC_COMM_WORLD,ts,ierr)
109:       call TSSetProblemType(ts,TS_NONLINEAR,ierr)

111: !
112: !     Tell the timestepper context where to compute solutions
113: !

115:       call TSSetSolution(ts,x,ierr)

117: !
118: !    Provide the call-back for the nonlinear function we are
119: !    evaluating. Thus whenever the timestepping routines need the
120: !    function they will call this routine. Note the final argument
121: !    is the application context used by the call-back functions.
122: !

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

126: !
127: !     Set the Jacobian matrix and the function used to compute
128: !     Jacobians.
129: !

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

133: !
134: !       For the initial guess for the problem
135: !

137:       call FormInitialGuess(x,user,ierr)

139: !
140: !       This indicates that we are using pseudo timestepping to
141: !     find a steady state solution to the nonlinear problem.
142: !

144:       call TSSetType(ts,TSPSEUDO,ierr)

146: !
147: !       Set the initial time to start at (this is arbitrary for
148: !     steady state problems and the initial timestep given above
149: !

151:       call TSSetTimeStep(ts,dt,ierr)

153: !
154: !      Set a large number of timesteps and final duration time
155: !     to insure convergence to steady state.
156: !
157:       i1000 = 1000
158:       tmax  = 1.e12
159:       call TSSetMaxSteps(ts,i1000,ierr)
160:       call TSSetMaxTime(ts,tmax,ierr)
161:       call TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr)

163: !
164: !      Set any additional options from the options database. This
165: !     includes all options for the nonlinear and linear solvers used
166: !     internally the timestepping routines.
167: !

169:       call TSSetFromOptions(ts,ierr)

171:       call TSSetUp(ts,ierr)

173: !
174: !      Perform the solve. This is where the timestepping takes place.
175: !
176:       call TSSolve(ts,x,ierr)
177:       call TSGetSolveTime(ts,ftime,ierr)
178:       call TSGetStepNumber(ts,its,ierr)
179:       call TSGetConvergedReason(ts,reason,ierr)

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

185: !
186: !     Free the data structures constructed above
187: !

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

196: !
197: !  --------------------  Form initial approximation -----------------
198: !
199:       subroutine FormInitialGuess(X,user,ierr)
200:       use petscts
201:       implicit none

203:       Vec              X
204:       PetscReal user(3)
205:       PetscInt  i,j,row,mx,my
206:       PetscErrorCode ierr
207:       PetscOffset      xidx
208:       PetscReal one,lambda
209:       PetscReal temp1,temp,hx,hy
210:       PetscScalar      xx(2)
211:       PetscInt          param,lmx,lmy
212:       parameter (param = 1,lmx = 2,lmy = 3)

214:       one = 1.0

216:       mx     = int(user(lmx))
217:       my     = int(user(lmy))
218:       lambda = user(param)

220:       hy    = one / (my-1)
221:       hx    = one / (mx-1)

223:       call VecGetArray(X,xx,xidx,ierr)
224:       temp1 = lambda/(lambda + one)
225:       do 10, j=1,my
226:         temp = min(j-1,my-j)*hy
227:         do 20 i=1,mx
228:           row = i + (j-1)*mx
229:           if (i .eq. 1 .or. j .eq. 1 .or.                               &
230:      &        i .eq. mx .or. j .eq. my) then
231:             xx(row+xidx) = 0.0
232:           else
233:             xx(row+xidx) =                                              &
234:      &        temp1*sqrt(min(min(i-1,mx-i)*hx,temp))
235:           endif
236:  20     continue
237:  10   continue
238:       call VecRestoreArray(X,xx,xidx,ierr)
239:       return
240:       end
241: !
242: !  --------------------  Evaluate Function F(x) ---------------------
243: !
244:       subroutine FormFunction(ts,t,X,F,user,ierr)
245:       use petscts
246:       implicit none

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

263:       two = 2.0

265:       mx     = int(user(lmx))
266:       my     = int(user(lmy))
267:       lambda = user(param)

269:       hx    = 1.0 / real(mx-1)
270:       hy    = 1.0 / real(my-1)
271:       sc    = hx*hy
272:       hxdhy = hx/hy
273:       hydhx = hy/hx

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

297:       call VecRestoreArrayRead(X,xx,xidx,ierr)
298:       call VecRestoreArray(F,ff,fidx,ierr)
299:       return
300:       end
301: !
302: !  --------------------  Evaluate Jacobian of F(x) --------------------
303: !
304:       subroutine FormJacobian(ts,ctime,X,JJ,B,user,ierr)
305:       use petscts
306:       implicit none

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

321:       PetscInt  param,lmx,lmy
322:       parameter (param = 1,lmx = 2,lmy = 3)

324:       i1 = 1
325:       i5 = 5
326:       jac = B
327:       two = 2.0
328:       one = 1.0

330:       mx     = int(user(lmx))
331:       my     = int(user(lmy))
332:       lambda = user(param)

334:       hx    = 1.0 / real(mx-1)
335:       hy    = 1.0 / real(my-1)
336:       sc    = hx*hy
337:       hxdhy = hx/hy
338:       hydhx = hy/hx

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

371: !/*TEST
372: !
373: !    test:
374: !      TODO: broken
375: !      args: -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_monitor_pseudo -ts_max_snes_failures 3 -ts_pseudo_frtol 1.e-5  -snes_stol 1e-5
376: !
377: !TEST*/