Actual source code: ex54f.F90
petsc-3.9.4 2018-09-11
1: !
2: ! Description: Solve Ax=b. A comes from an anisotropic 2D thermal problem with Q1 FEM on domain (-1,1)^2.
3: ! Material conductivity given by tensor:
4: !
5: ! D = | 1 0 |
6: ! | 0 epsilon |
7: !
8: ! rotated by angle 'theta' (-theta <90> in degrees) with anisotropic parameter 'epsilon' (-epsilon <0.0>).
9: ! Blob right hand side centered at C (-blob_center C(1),C(2) <0,0>)
10: ! Dirichlet BCs on y=-1 face.
11: !
12: ! -out_matlab will generate binary files for A,x,b and a ex54f.m file that reads them and plots them in matlab.
13: !
14: ! User can change anisotropic shape with function ex54_psi(). Negative theta will switch to a circular anisotropy.
15: !
16: !!/*T
17: ! Concepts: KSP^solving a system of linear equations
18: !! requires: !single
19: !T*/
22: ! -----------------------------------------------------------------------
23: program main
24: #include <petsc/finclude/petscksp.h>
25: use petscksp
26: implicit none
28: Vec xvec,bvec,uvec
29: Mat Amat
30: KSP ksp
31: PetscErrorCode ierr
32: PetscViewer viewer
33: PetscInt qj,qi,ne,M,Istart,Iend,geq,ix
34: PetscInt ki,kj,lint,nel,ll,j1,i1,ndf,f4
35: PetscInt f2,f9,f6,one
36: PetscInt :: idx(4)
37: PetscBool flg,out_matlab
38: PetscMPIInt size,rank
39: PetscScalar::ss(4,4),val
40: PetscReal::shp(3,9),sg(3,9)
41: PetscReal::thk,a1,a2
42: PetscReal, external :: ex54_psi
43: PetscReal::theta,eps,h,x,y,xsj
44: PetscReal::coord(2,4),dd(2,2),ev(3),blb(2)
46: common /ex54_theta/ theta
47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: ! Beginning of program
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
51: if (ierr .ne. 0) then
52: print*,'Unable to initialize PETSc'
53: stop
54: endif
55: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
56: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
57: one = 1
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: ! set parameters
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: f4 = 4
62: f2 = 2
63: f9 = 9
64: f6 = 6
65: ne = 9
66: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-ne',ne,flg,ierr)
67: h = 2.0/real(ne)
68: M = (ne+1)*(ne+1)
69: theta = 90.0
70: ! theta is input in degrees
71: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-theta',theta,flg,ierr)
72: theta = theta / 57.2957795
73: eps = 1.0
74: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-epsilon',eps,flg,ierr)
75: ki = 2
76: call PetscOptionsGetRealArray(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-blob_center',blb,ki,flg,ierr)
77: if ( .not. flg ) then
78: blb(1) = 0.0
79: blb(2) = 0.0
80: else if ( ki .ne. 2 ) then
81: print *, 'error: ', ki,' arguments read for -blob_center. Needs to be two.'
82: endif
83: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-out_matlab',out_matlab,flg,ierr)
84: if (.not.flg) out_matlab = PETSC_FALSE;
86: ev(1) = 1.0
87: ev(2) = eps*ev(1)
88: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: ! Compute the matrix and right-hand-side vector that define
90: ! the linear system, Ax = b.
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: ! Create matrix. When using MatCreate(), the matrix format can
93: ! be specified at runtime.
94: call MatCreate(PETSC_COMM_WORLD,Amat,ierr)
95: call MatSetSizes( Amat,PETSC_DECIDE, PETSC_DECIDE, M, M, ierr )
96: if ( size == 1 ) then
97: call MatSetType( Amat, MATAIJ, ierr )
98: else
99: call MatSetType( Amat, MATMPIAIJ, ierr )
100: endif
101: call MatMPIAIJSetPreallocation(Amat,f9,PETSC_NULL_INTEGER,f6,PETSC_NULL_INTEGER, ierr)
102: call MatSetFromOptions( Amat, ierr )
103: call MatSetUp( Amat, ierr )
104: call MatGetOwnershipRange( Amat, Istart, Iend, ierr )
105: ! Create vectors. Note that we form 1 vector from scratch and
106: ! then duplicate as needed.
107: xvec = tVec(0)
108: call MatCreateVecs( Amat, PETSC_NULL_VEC, xvec, ierr )
109: call VecSetFromOptions( xvec, ierr )
110: call VecDuplicate( xvec, bvec, ierr )
111: call VecDuplicate( xvec, uvec, ierr )
112: ! Assemble matrix.
113: ! - Note that MatSetValues() uses 0-based row and column numbers
114: ! in Fortran as well as in C (as set here in the array "col").
115: thk = 1.0 ! thickness
116: nel = 4 ! nodes per element (quad)
117: ndf = 1
118: call int2d(f2,sg)
119: lint = 4
120: ix = 0
121: do geq=Istart,Iend-1,1
122: qj = geq/(ne+1); qi = mod(geq,(ne+1))
123: x = h*qi - 1.0; y = h*qj - 1.0 ! lower left corner (-1,-1)
124: if ( qi < ne .and. qj < ne ) then
125: coord(1,1) = x; coord(2,1) = y
126: coord(1,2) = x+h; coord(2,2) = y
127: coord(1,3) = x+h; coord(2,3) = y+h
128: coord(1,4) = x; coord(2,4) = y+h
129: ! form stiff
130: ss = 0.0
131: do ll = 1,lint
132: call shp2dquad(sg(1,ll),sg(2,ll),coord,shp,xsj,f2)
133: xsj = xsj*sg(3,ll)*thk
134: call thfx2d(ev,coord,shp,dd,f2,f2,f4,ex54_psi)
135: j1 = 1
136: do kj = 1,nel
137: a1 = (dd(1,1)*shp(1,kj) + dd(1,2)*shp(2,kj))*xsj
138: a2 = (dd(2,1)*shp(1,kj) + dd(2,2)*shp(2,kj))*xsj
139: ! Compute residual
140: ! p(j1) = p(j1) - a1*gradt(1) - a2*gradt(2)
141: ! Compute tangent
142: i1 = 1
143: do ki = 1,nel
144: ss(i1,j1) = ss(i1,j1) + a1*shp(1,ki) + a2*shp(2,ki)
145: i1 = i1 + ndf
146: end do
147: j1 = j1 + ndf
148: end do
149: enddo
151: idx(1) = geq; idx(2) = geq+1; idx(3) = geq+(ne+1)+1
152: idx(4) = geq+(ne+1)
153: if ( qj > 0 ) then
154: call MatSetValues(Amat,f4,idx,f4,idx,ss,ADD_VALUES,ierr)
155: else ! a BC
156: do ki=1,4,1
157: do kj=1,4,1
158: if (ki<3 .or. kj<3 ) then
159: if ( ki==kj ) then
160: ss(ki,kj) = .1*ss(ki,kj)
161: else
162: ss(ki,kj) = 0.0
163: endif
164: endif
165: enddo
166: enddo
167: call MatSetValues(Amat,f4,idx,f4,idx,ss,ADD_VALUES,ierr)
168: endif ! BC
169: endif ! add element
170: if ( qj > 0 ) then ! set rhs
171: val = h*h*exp(-100*((x+h/2)-blb(1))**2)*exp(-100*((y+h/2)-blb(2))**2)
172: call VecSetValues(bvec,one,geq,val,INSERT_VALUES,ierr)
173: endif
174: enddo
175: call MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY,ierr)
176: call MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY,ierr)
177: call VecAssemblyBegin(bvec,ierr)
178: call VecAssemblyEnd(bvec,ierr)
180: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: ! Create the linear solver and set various options
182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: ! Create linear solver context
186: call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
188: ! Set operators. Here the matrix that defines the linear system
189: ! also serves as the preconditioning matrix.
191: call KSPSetOperators(ksp,Amat,Amat,ierr)
193: ! Set runtime options, e.g.,
194: ! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
195: ! These options will override those specified above as long as
196: ! KSPSetFromOptions() is called _after_ any other customization
197: ! routines.
199: call KSPSetFromOptions(ksp,ierr)
201: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: ! Solve the linear system
203: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: call KSPSolve(ksp,bvec,xvec,ierr)
206: CHKERRA(ierr)
209: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210: ! output
211: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212: if ( out_matlab ) then
213: call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Amat',FILE_MODE_WRITE,viewer,ierr)
214: call MatView(Amat,viewer,ierr)
215: call PetscViewerDestroy(viewer,ierr)
217: call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Bvec',FILE_MODE_WRITE,viewer,ierr)
218: call VecView(bvec,viewer,ierr)
219: call PetscViewerDestroy(viewer,ierr)
221: call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Xvec',FILE_MODE_WRITE,viewer,ierr)
222: call VecView(xvec,viewer,ierr)
223: call PetscViewerDestroy(viewer,ierr)
225: call MatMult(Amat,xvec,uvec,ierr)
226: val = -1.0
227: call VecAXPY(uvec,val,bvec,ierr)
228: call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Rvec',FILE_MODE_WRITE,viewer,ierr)
229: call VecView(uvec,viewer,ierr)
230: call PetscViewerDestroy(viewer,ierr)
232: if ( rank == 0 ) then
233: open(1,file='ex54f.m', FORM='formatted')
234: write (1,*) 'A = PetscBinaryRead(''Amat'');'
235: write (1,*) '[m n] = size(A);'
236: write (1,*) 'mm = sqrt(m);'
237: write (1,*) 'b = PetscBinaryRead(''Bvec'');'
238: write (1,*) 'x = PetscBinaryRead(''Xvec'');'
239: write (1,*) 'r = PetscBinaryRead(''Rvec'');'
240: write (1,*) 'bb = reshape(b,mm,mm);'
241: write (1,*) 'xx = reshape(x,mm,mm);'
242: write (1,*) 'rr = reshape(r,mm,mm);'
243: ! write (1,*) 'imagesc(bb')'
244: ! write (1,*) 'title('RHS'),'
245: write (1,*) 'figure,'
246: write (1,*) 'imagesc(xx'')'
247: write (1,2002) eps,theta*57.2957795
248: write (1,*) 'figure,'
249: write (1,*) 'imagesc(rr'')'
250: write (1,*) 'title(''Residual''),'
251: close(1)
252: endif
253: endif
254: 2002 format('title(''Solution: esp='',d9.3,'', theta='',g8.3,''),')
255: ! Free work space. All PETSc objects should be destroyed when they
256: ! are no longer needed.
258: call VecDestroy(xvec,ierr)
259: call VecDestroy(bvec,ierr)
260: call VecDestroy(uvec,ierr)
261: call MatDestroy(Amat,ierr)
262: call KSPDestroy(ksp,ierr)
263: call PetscFinalize(ierr)
265: end
267: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268: ! thfx2d - compute material tensor
269: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: ! Compute thermal gradient and flux
272: subroutine thfx2d(ev,xl,shp,dd,ndm,ndf,nel,dir)
273: implicit none
275: PetscInt ndm,ndf,nel,i
276: PetscReal ev(2),xl(ndm,nel),shp(3,*),dir
277: PetscReal xx,yy,psi,cs,sn,c2,s2,dd(2,2)
279: xx = 0.0
280: yy = 0.0
281: do i = 1,nel
282: xx = xx + shp(3,i)*xl(1,i)
283: yy = yy + shp(3,i)*xl(2,i)
284: end do
285: psi = dir(xx,yy)
286: ! Compute thermal flux
287: cs = cos(psi)
288: sn = sin(psi)
289: c2 = cs*cs
290: s2 = sn*sn
291: cs = cs*sn
293: dd(1,1) = c2*ev(1) + s2*ev(2)
294: dd(2,2) = s2*ev(1) + c2*ev(2)
295: dd(1,2) = cs*(ev(1) - ev(2))
296: dd(2,1) = dd(1,2)
298: ! flux(1) = -dd(1,1)*gradt(1) - dd(1,2)*gradt(2)
299: ! flux(2) = -dd(2,1)*gradt(1) - dd(2,2)*gradt(2)
301: end
303: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
304: ! shp2dquad - shape functions - compute derivatives w/r natural coords.
305: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
306: subroutine shp2dquad(s,t,xl,shp,xsj,ndm)
307: !-----[--.----+----.----+----.-----------------------------------------]
308: ! Purpose: Shape function routine for 4-node isoparametric quads
309: !
310: ! Inputs:
311: ! s,t - Natural coordinates of point
312: ! xl(ndm,*) - Nodal coordinates for element
313: ! ndm - Spatial dimension of mesh
315: ! Outputs:
316: ! shp(3,*) - Shape functions and derivatives at point
317: ! shp(1,i) = dN_i/dx or dN_i/dxi_1
318: ! shp(2,i) = dN_i/dy or dN_i/dxi_2
319: ! shp(3,i) = N_i
320: ! xsj - Jacobian determinant at point
321: !-----[--.----+----.----+----.-----------------------------------------]
322: implicit none
323: PetscInt ndm
324: PetscReal xo,xs,xt, yo,ys,yt, xsm,xsp,xtm
325: PetscReal xtp, ysm,ysp,ytm,ytp
326: PetscReal s,t, xsj,xsj1, sh,th,sp,tp,sm
327: PetscReal tm, xl(ndm,4),shp(3,4)
329: ! Set up interpolations
331: sh = 0.5*s
332: th = 0.5*t
333: sp = 0.5 + sh
334: tp = 0.5 + th
335: sm = 0.5 - sh
336: tm = 0.5 - th
337: shp(3,1) = sm*tm
338: shp(3,2) = sp*tm
339: shp(3,3) = sp*tp
340: shp(3,4) = sm*tp
342: ! Set up natural coordinate functions (times 4)
344: xo = xl(1,1)-xl(1,2)+xl(1,3)-xl(1,4)
345: xs = -xl(1,1)+xl(1,2)+xl(1,3)-xl(1,4) + xo*t
346: xt = -xl(1,1)-xl(1,2)+xl(1,3)+xl(1,4) + xo*s
347: yo = xl(2,1)-xl(2,2)+xl(2,3)-xl(2,4)
348: ys = -xl(2,1)+xl(2,2)+xl(2,3)-xl(2,4) + yo*t
349: yt = -xl(2,1)-xl(2,2)+xl(2,3)+xl(2,4) + yo*s
351: ! Compute jacobian (times 16)
353: xsj1 = xs*yt - xt*ys
355: ! Divide jacobian by 16 (multiply by .0625)
357: xsj = 0.0625*xsj1
358: if (xsj1.eq.0.0) then
359: xsj1 = 1.0
360: else
361: xsj1 = 1.0/xsj1
362: endif
364: ! Divide functions by jacobian
366: xs = (xs+xs)*xsj1
367: xt = (xt+xt)*xsj1
368: ys = (ys+ys)*xsj1
369: yt = (yt+yt)*xsj1
371: ! Multiply by interpolations
373: ytm = yt*tm
374: ysm = ys*sm
375: ytp = yt*tp
376: ysp = ys*sp
377: xtm = xt*tm
378: xsm = xs*sm
379: xtp = xt*tp
380: xsp = xs*sp
382: ! Compute shape functions
384: shp(1,1) = - ytm+ysm
385: shp(1,2) = ytm+ysp
386: shp(1,3) = ytp-ysp
387: shp(1,4) = - ytp-ysm
388: shp(2,1) = xtm-xsm
389: shp(2,2) = - xtm-xsp
390: shp(2,3) = - xtp+xsp
391: shp(2,4) = xtp+xsm
393: end
395: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396: ! int2d
397: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
398: subroutine int2d(l,sg)
399: !-----[--.----+----.----+----.-----------------------------------------]
400: ! Purpose: Form Gauss points and weights for two dimensions
402: ! Inputs:
403: ! l - Number of points/direction
405: ! Outputs:
406: ! lint - Total number of points
407: ! sg(3,*) - Array of points and weights
408: !-----[--.----+----.----+----.-----------------------------------------]
409: implicit none
410: PetscInt l,i,lint,lr(9),lz(9)
411: PetscReal g,third,sg(3,*)
412: data lr/-1,1,1,-1,0,1,0,-1,0/,lz/-1,-1,1,1,-1,0,1,0,0/
413: data third / 0.3333333333333333 /
415: ! Set number of total points
417: lint = l*l
419: ! 2x2 integration
420: g = sqrt(third)
421: do i = 1,4
422: sg(1,i) = g*lr(i)
423: sg(2,i) = g*lz(i)
424: sg(3,i) = 1.0
425: end do
427: end
429: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
430: ! ex54_psi - anusotropic material direction
431: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432: PetscReal function ex54_psi(x,y)
433: implicit none
434: PetscReal x,y,theta
435: common /ex54_theta/ theta
436: ex54_psi = theta
437: if ( theta < 0. ) then ! circular
438: if (y==0) then
439: ex54_psi = 2.0*atan(1.0)
440: else
441: ex54_psi = atan(-x/y)
442: endif
443: endif
444: end
446: !
447: !/*TEST
448: !
449: ! build:
450: ! requires: !libpgf90
451: !
452: ! test:
453: ! nsize: 4
454: ! args: -ne 39 -theta 30.0 -epsilon 1.e-1 -blob_center 0.,0. -ksp_type cg -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05 -mat_coarsen_type hem -pc_gamg_square_graph 0 -ksp_monitor_short -mg_levels_esteig_ksp_type cg
455: ! requires: !single
456: !
457: !TEST*/