Actual source code: ex1f.F90

  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:       PetscCallA(PetscInitialize(ierr))
 56:       user(lmx)        = 4
 57:       user(lmy)        = 4
 58:       user(param)      = 6.0

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

 79: !
 80: !      Create vectors to hold the solution and function value
 81: !
 82:       PetscCallA(VecCreateSeq(PETSC_COMM_SELF,N,x,ierr))
 83:       PetscCallA(VecDuplicate(x,r,ierr))

 85: !
 86: !    Create matrix to hold Jacobian. Preallocate 5 nonzeros per row
 87: !    in the sparse matrix. Note that this is not the optimal strategy see
 88: !    the Performance chapter of the users manual for information on
 89: !    preallocating memory in sparse matrices.
 90: !
 91:       i5 = 5
 92:       PetscCallA(MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,i5,PETSC_NULL_INTEGER,J,ierr))

 94: !
 95: !     Create timestepper context
 96: !

 98:       PetscCallA(TSCreate(PETSC_COMM_WORLD,ts,ierr))
 99:       PetscCallA(TSSetProblemType(ts,TS_NONLINEAR,ierr))

101: !
102: !     Tell the timestepper context where to compute solutions
103: !

105:       PetscCallA(TSSetSolution(ts,x,ierr))

107: !
108: !    Provide the call-back for the nonlinear function we are
109: !    evaluating. Thus whenever the timestepping routines need the
110: !    function they will call this routine. Note the final argument
111: !    is the application context used by the call-back functions.
112: !

114:       PetscCallA(TSSetRHSFunction(ts,PETSC_NULL_VEC,FormFunction,user,ierr))

116: !
117: !     Set the Jacobian matrix and the function used to compute
118: !     Jacobians.
119: !

121:       PetscCallA(TSSetRHSJacobian(ts,J,J,FormJacobian,user,ierr))

123: !
124: !       For the initial guess for the problem
125: !

127:       PetscCallA(FormInitialGuess(x,user,ierr))

129: !
130: !       This indicates that we are using pseudo timestepping to
131: !     find a steady state solution to the nonlinear problem.
132: !

134:       PetscCallA(TSSetType(ts,TSPSEUDO,ierr))

136: !
137: !       Set the initial time to start at (this is arbitrary for
138: !     steady state problems and the initial timestep given above
139: !

141:       PetscCallA(TSSetTimeStep(ts,dt,ierr))

143: !
144: !      Set a large number of timesteps and final duration time
145: !     to insure convergence to steady state.
146: !
147:       i1000 = 1000
148:       tmax  = 1.e12
149:       PetscCallA(TSSetMaxSteps(ts,i1000,ierr))
150:       PetscCallA(TSSetMaxTime(ts,tmax,ierr))
151:       PetscCallA(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER,ierr))

153: !
154: !      Set any additional options from the options database. This
155: !     includes all options for the nonlinear and linear solvers used
156: !     internally the timestepping routines.
157: !

159:       PetscCallA(TSSetFromOptions(ts,ierr))

161:       PetscCallA(TSSetUp(ts,ierr))

163: !
164: !      Perform the solve. This is where the timestepping takes place.
165: !
166:       PetscCallA(TSSolve(ts,x,ierr))
167:       PetscCallA(TSGetSolveTime(ts,ftime,ierr))
168:       PetscCallA(TSGetStepNumber(ts,its,ierr))
169:       PetscCallA(TSGetConvergedReason(ts,reason,ierr))

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

174: !
175: !     Free the data structures constructed above
176: !

178:       PetscCallA(VecDestroy(x,ierr))
179:       PetscCallA(VecDestroy(r,ierr))
180:       PetscCallA(MatDestroy(J,ierr))
181:       PetscCallA(TSDestroy(ts,ierr))
182:       PetscCallA(PetscFinalize(ierr))
183:       end

185: !
186: !  --------------------  Form initial approximation -----------------
187: !
188:       subroutine FormInitialGuess(X,user,ierr)
189:       use petscts
190:       implicit none

192:       Vec              X
193:       PetscReal user(3)
194:       PetscInt  i,j,row,mx,my
195:       PetscErrorCode ierr
196:       PetscOffset      xidx
197:       PetscReal one,lambda
198:       PetscReal temp1,temp,hx,hy
199:       PetscScalar      xx(2)
200:       PetscInt          param,lmx,lmy
201:       parameter (param = 1,lmx = 2,lmy = 3)

203:       one = 1.0

205:       mx     = int(user(lmx))
206:       my     = int(user(lmy))
207:       lambda = user(param)

209:       hy    = one / (my-1)
210:       hx    = one / (mx-1)

212:       PetscCall(VecGetArray(X,xx,xidx,ierr))
213:       temp1 = lambda/(lambda + one)
214:       do 10, j=1,my
215:         temp = min(j-1,my-j)*hy
216:         do 20 i=1,mx
217:           row = i + (j-1)*mx
218:           if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
219:             xx(row+xidx) = 0.0
220:           else
221:             xx(row+xidx) = temp1*sqrt(min(min(i-1,mx-i)*hx,temp))
222:           endif
223:  20     continue
224:  10   continue
225:       PetscCall(VecRestoreArray(X,xx,xidx,ierr))
226:       return
227:       end
228: !
229: !  --------------------  Evaluate Function F(x) ---------------------
230: !
231:       subroutine FormFunction(ts,t,X,F,user,ierr)
232:       use petscts
233:       implicit none

235:       TS       ts
236:       PetscReal  t
237:       Vec               X,F
238:       PetscReal  user(3)
239:       PetscErrorCode     ierr
240:       PetscInt         i,j,row,mx,my
241:       PetscOffset       xidx,fidx
242:       PetscReal  two,lambda
243:       PetscReal  hx,hy,hxdhy,hydhx
244:       PetscScalar  ut,ub,ul,ur,u
245:       PetscScalar  uxx,uyy,sc
246:       PetscScalar  xx(2),ff(2)
247:       PetscInt     param,lmx,lmy
248:       parameter (param = 1,lmx = 2,lmy = 3)

250:       two = 2.0

252:       mx     = int(user(lmx))
253:       my     = int(user(lmy))
254:       lambda = user(param)

256:       hx    = 1.0 / real(mx-1)
257:       hy    = 1.0 / real(my-1)
258:       sc    = hx*hy
259:       hxdhy = hx/hy
260:       hydhx = hy/hx

262:       PetscCall(VecGetArrayRead(X,xx,xidx,ierr))
263:       PetscCall(VecGetArray(F,ff,fidx,ierr))
264:       do 10 j=1,my
265:         do 20 i=1,mx
266:           row = i + (j-1)*mx
267:           if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
268:             ff(row+fidx) = xx(row+xidx)
269:           else
270:             u            = xx(row + xidx)
271:             ub           = xx(row - mx + xidx)
272:             ul           = xx(row - 1 + xidx)
273:             ut           = xx(row + mx + xidx)
274:             ur           = xx(row + 1 + xidx)
275:             uxx          = (-ur + two*u - ul)*hydhx
276:             uyy          = (-ut + two*u - ub)*hxdhy
277:             ff(row+fidx) = -uxx - uyy + sc*lambda*exp(u)
278:             u =  -uxx - uyy + sc*lambda*exp(u)
279:          endif
280:  20   continue
281:  10   continue

283:       PetscCall(VecRestoreArrayRead(X,xx,xidx,ierr))
284:       PetscCall(VecRestoreArray(F,ff,fidx,ierr))
285:       return
286:       end
287: !
288: !  --------------------  Evaluate Jacobian of F(x) --------------------
289: !
290:       subroutine FormJacobian(ts,ctime,X,JJ,B,user,ierr)
291:       use petscts
292:       implicit none

294:       TS               ts
295:       Vec              X
296:       Mat              JJ,B
297:       PetscReal user(3),ctime
298:       PetscErrorCode   ierr
299:       Mat              jac
300:       PetscOffset xidx
301:       PetscInt    i,j,row(1),mx,my
302:       PetscInt    col(5),i1,i5
303:       PetscScalar two,one,lambda
304:       PetscScalar v(5),sc,xx(2)
305:       PetscReal hx,hy,hxdhy,hydhx

307:       PetscInt  param,lmx,lmy
308:       parameter (param = 1,lmx = 2,lmy = 3)

310:       i1 = 1
311:       i5 = 5
312:       jac = B
313:       two = 2.0
314:       one = 1.0

316:       mx     = int(user(lmx))
317:       my     = int(user(lmy))
318:       lambda = user(param)

320:       hx    = 1.0 / real(mx-1)
321:       hy    = 1.0 / real(my-1)
322:       sc    = hx*hy
323:       hxdhy = hx/hy
324:       hydhx = hy/hx

326:       PetscCall(VecGetArrayRead(X,xx,xidx,ierr))
327:       do 10 j=1,my
328:         do 20 i=1,mx
329: !
330: !      When inserting into PETSc matrices, indices start at 0
331: !
332:           row(1) = i - 1 + (j-1)*mx
333:           if (i .eq. 1 .or. j .eq. 1 .or. i .eq. mx .or. j .eq. my) then
334:             PetscCall(MatSetValues(jac,i1,row,i1,row,one,INSERT_VALUES,ierr))
335:           else
336:             v(1)   = hxdhy
337:             col(1) = row(1) - mx
338:             v(2)   = hydhx
339:             col(2) = row(1) - 1
340:             v(3)   = -two*(hydhx+hxdhy)+sc*lambda*exp(xx(row(1)+1+xidx))
341:             col(3) = row(1)
342:             v(4)   = hydhx
343:             col(4) = row(1) + 1
344:             v(5)   = hxdhy
345:             col(5) = row(1) + mx
346:             PetscCall(MatSetValues(jac,i1,row,i5,col,v,INSERT_VALUES,ierr))
347:           endif
348:  20     continue
349:  10   continue
350:       PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
351:       PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
352:       PetscCall(VecRestoreArray(X,xx,xidx,ierr))
353:       return
354:       end

356: !/*TEST
357: !
358: !    test:
359: !      TODO: broken
360: !      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
361: !
362: !TEST*/