Actual source code: ex1f.F
petsc-3.14.6 2021-03-30
1: !
2: ! Solves the time dependent Bratu problem using pseudo-timestepping
3: !
4: ! Concepts: TS^pseudo-timestepping
5: ! Concepts: TS^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/tutorials/ex4.c[ex4f.F] and
19: ! snes/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
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: call TSSetTimeStep(ts,dt,ierr)
158: !
159: ! Set a large number of timesteps and final duration time
160: ! to insure convergence to steady state.
161: !
162: i1000 = 1000
163: tmax = 1.e12
164: call TSSetMaxSteps(ts,i1000,ierr)
165: call TSSetMaxTime(ts,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 TSGetStepNumber(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 ',1pe9.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
376: !/*TEST
377: !
378: ! test:
379: ! TODO: broken
380: ! 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
381: !
382: !TEST*/