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*/