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