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