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