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