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: !
 16: !!/*T
 17: !   Concepts: KSP^solving a system of linear equations
 18: !T*/

 20: ! -----------------------------------------------------------------------
 21:       program main
 22: #include <petsc/finclude/petscksp.h>
 23:       use petscksp
 24:       implicit none

 26:       Vec              xvec,bvec,uvec
 27:       Mat              Amat
 28:       KSP              ksp
 29:       PetscErrorCode   ierr
 30:       PetscViewer viewer
 31:       PetscInt qj,qi,ne,M,Istart,Iend,geq
 32:       PetscInt ki,kj,nel,ll,j1,i1,ndf,f4
 33:       PetscInt f2,f9,f6,one
 34:       PetscInt :: idx(4)
 35:       PetscBool  flg,out_matlab
 36:       PetscMPIInt size,rank
 37:       PetscScalar::ss(4,4),val
 38:       PetscReal::shp(3,9),sg(3,9)
 39:       PetscReal::thk,a1,a2
 40:       PetscReal, external :: ex54_psi
 41:       PetscReal::theta,eps,h,x,y,xsj
 42:       PetscReal::coord(2,4),dd(2,2),ev(3),blb(2)

 44:       common /ex54_theta/ theta
 45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 46: !                 Beginning of program
 47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 49:       if (ierr .ne. 0) then
 50:         print*,'Unable to initialize PETSc'
 51:         stop
 52:       endif
 53:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 54:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 55:       one = 1
 56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 57: !                 set parameters
 58: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 59:       f4 = 4
 60:       f2 = 2
 61:       f9 = 9
 62:       f6 = 6
 63:       ne = 9
 64:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-ne',ne,flg,ierr)
 65:       h = 2.0/real(ne)
 66:       M = (ne+1)*(ne+1)
 67:       theta = 90.0
 68: !     theta is input in degrees
 69:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-theta',theta,flg,ierr)
 70:       theta = theta / 57.2957795
 71:       eps = 1.0
 72:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-epsilon',eps,flg,ierr)
 73:       ki = 2
 74:       call PetscOptionsGetRealArray(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-blob_center',blb,ki,flg,ierr)
 75:       if (.not. flg) then
 76:          blb(1) = 0.0
 77:          blb(2) = 0.0
 78:       else if (ki .ne. 2) then
 79:          print *, 'error: ', ki,' arguments read for -blob_center.  Needs to be two.'
 80:       endif
 81:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-out_matlab',out_matlab,flg,ierr)
 82:       if (.not.flg) out_matlab = PETSC_FALSE;

 84:       ev(1) = 1.0
 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:       call MatSetType( Amat, MATAIJ, ierr)
 95:       call MatSetOption(Amat,MAT_SPD,PETSC_TRUE,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:       call MatCreateVecs( Amat, PETSC_NULL_VEC, xvec, ierr)
108:       call VecSetFromOptions( xvec, ierr)
109:       call VecDuplicate( xvec, bvec, ierr)
110:       call VecDuplicate( xvec, uvec, ierr)
111: !  Assemble matrix.
112: !   - Note that MatSetValues() uses 0-based row and column numbers
113: !     in Fortran as well as in C (as set here in the array "col").
114:       thk = 1.0              ! thickness
115:       nel = 4                   ! nodes per element (quad)
116:       ndf = 1
117:       call int2d(f2,sg)
118:       do geq=Istart,Iend-1,1
119:          qj = geq/(ne+1); qi = mod(geq,(ne+1))
120:          x = h*qi - 1.0; y = h*qj - 1.0 ! 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.0
128:             do ll = 1,4
129:                call shp2dquad(sg(1,ll),sg(2,ll),coord,shp,xsj,f2)
130:                xsj = xsj*sg(3,ll)*thk
131:                call thfx2d(ev,coord,shp,dd,f2,f2,f4,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: !     Compute residual
137: !                  p(j1) = p(j1) - a1*gradt(1) - a2*gradt(2)
138: !     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,f4,idx,f4,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) = .1*ss(ki,kj)
158:                         else
159:                            ss(ki,kj) = 0.0
160:                         endif
161:                      endif
162:                   enddo
163:                enddo
164:                call MatSetValues(Amat,f4,idx,f4,idx,ss,ADD_VALUES,ierr)
165:             endif               ! BC
166:          endif                  ! add element
167:          if (qj > 0) then      ! set rhs
168:             val = h*h*exp(-100*((x+h/2)-blb(1))**2)*exp(-100*((y+h/2)-blb(2))**2)
169:             call VecSetValues(bvec,one,geq,val,INSERT_VALUES,ierr)
170:          endif
171:       enddo
172:       call MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY,ierr)
173:       call MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY,ierr)
174:       call VecAssemblyBegin(bvec,ierr)
175:       call VecAssemblyEnd(bvec,ierr)

177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: !          Create the linear solver and set various options
179: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

181: !  Create linear solver context

183:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

185: !  Set operators. Here the matrix that defines the linear system
186: !  also serves as the preconditioning matrix.

188:       call KSPSetOperators(ksp,Amat,Amat,ierr)

190: !  Set runtime options, e.g.,
191: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
192: !  These options will override those specified above as long as
193: !  KSPSetFromOptions() is called _after_ any other customization
194: !  routines.

196:       call KSPSetFromOptions(ksp,ierr)

198: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: !                      Solve the linear system
200: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

202:       call KSPSolve(ksp,bvec,xvec,ierr)
203:       CHKERRA(ierr)

205: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206: !                      output
207: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208:       if (out_matlab) then
209:          call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Amat',FILE_MODE_WRITE,viewer,ierr)
210:          call MatView(Amat,viewer,ierr)
211:          call PetscViewerDestroy(viewer,ierr)

213:          call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Bvec',FILE_MODE_WRITE,viewer,ierr)
214:          call VecView(bvec,viewer,ierr)
215:          call PetscViewerDestroy(viewer,ierr)

217:          call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Xvec',FILE_MODE_WRITE,viewer,ierr)
218:          call VecView(xvec,viewer,ierr)
219:          call PetscViewerDestroy(viewer,ierr)

221:          call MatMult(Amat,xvec,uvec,ierr)
222:          val = -1.0
223:          call VecAXPY(uvec,val,bvec,ierr)
224:          call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Rvec',FILE_MODE_WRITE,viewer,ierr)
225:          call VecView(uvec,viewer,ierr)
226:          call PetscViewerDestroy(viewer,ierr)

228:          if (rank == 0) then
229:             open(1,file='ex54f.m', FORM='formatted')
230:             write (1,*) 'A = PetscBinaryRead(''Amat'');'
231:             write (1,*) '[m n] = size(A);'
232:             write (1,*) 'mm = sqrt(m);'
233:             write (1,*) 'b = PetscBinaryRead(''Bvec'');'
234:             write (1,*) 'x = PetscBinaryRead(''Xvec'');'
235:             write (1,*) 'r = PetscBinaryRead(''Rvec'');'
236:             write (1,*) 'bb = reshape(b,mm,mm);'
237:             write (1,*) 'xx = reshape(x,mm,mm);'
238:             write (1,*) 'rr = reshape(r,mm,mm);'
239: !            write (1,*) 'imagesc(bb')'
240: !            write (1,*) 'title('RHS'),'
241:             write (1,*) 'figure,'
242:             write (1,*) 'imagesc(xx'')'
243:             write (1,2002) eps,theta*57.2957795
244:             write (1,*) 'figure,'
245:             write (1,*) 'imagesc(rr'')'
246:             write (1,*) 'title(''Residual''),'
247:             close(1)
248:          endif
249:       endif
250:  2002 format('title(''Solution: esp='',d9.3,'', theta='',g8.3,''),')
251: !  Free work space.  All PETSc objects should be destroyed when they
252: !  are no longer needed.

254:       call VecDestroy(xvec,ierr)
255:       call VecDestroy(bvec,ierr)
256:       call VecDestroy(uvec,ierr)
257:       call MatDestroy(Amat,ierr)
258:       call KSPDestroy(ksp,ierr)
259:       call PetscFinalize(ierr)

261:       end

263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: !     thfx2d - compute material tensor
265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266: !     Compute thermal gradient and flux

268:       subroutine thfx2d(ev,xl,shp,dd,ndm,ndf,nel,dir)
269:       implicit  none

271:       PetscInt   ndm,ndf,nel,i
272:       PetscReal ev(2),xl(ndm,nel),shp(3,*),dir
273:       PetscReal xx,yy,psi,cs,sn,c2,s2,dd(2,2)

275:       xx       = 0.0
276:       yy       = 0.0
277:       do i = 1,nel
278:         xx       = xx       + shp(3,i)*xl(1,i)
279:         yy       = yy       + shp(3,i)*xl(2,i)
280:       end do
281:       psi = dir(xx,yy)
282: !     Compute thermal flux
283:       cs  = cos(psi)
284:       sn  = sin(psi)
285:       c2  = cs*cs
286:       s2  = sn*sn
287:       cs  = cs*sn

289:       dd(1,1) = c2*ev(1) + s2*ev(2)
290:       dd(2,2) = s2*ev(1) + c2*ev(2)
291:       dd(1,2) = cs*(ev(1) - ev(2))
292:       dd(2,1) = dd(1,2)

294: !      flux(1) = -dd(1,1)*gradt(1) - dd(1,2)*gradt(2)
295: !      flux(2) = -dd(2,1)*gradt(1) - dd(2,2)*gradt(2)

297:       end

299: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
300: !     shp2dquad - shape functions - compute derivatives w/r natural coords.
301: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
302:        subroutine shp2dquad(s,t,xl,shp,xsj,ndm)
303: !-----[--.----+----.----+----.-----------------------------------------]
304: !      Purpose: Shape function routine for 4-node isoparametric quads
305: !
306: !      Inputs:
307: !         s,t       - Natural coordinates of point
308: !         xl(ndm,*) - Nodal coordinates for element
309: !         ndm       - Spatial dimension of mesh

311: !      Outputs:
312: !         shp(3,*)  - Shape functions and derivatives at point
313: !                     shp(1,i) = dN_i/dx  or dN_i/dxi_1
314: !                     shp(2,i) = dN_i/dy  or dN_i/dxi_2
315: !                     shp(3,i) = N_i
316: !         xsj       - Jacobian determinant at point
317: !-----[--.----+----.----+----.-----------------------------------------]
318:       implicit  none
319:       PetscInt  ndm
320:       PetscReal xo,xs,xt, yo,ys,yt, xsm,xsp,xtm
321:       PetscReal xtp, ysm,ysp,ytm,ytp
322:       PetscReal s,t, xsj,xsj1, sh,th,sp,tp,sm
323:       PetscReal tm, xl(ndm,4),shp(3,4)

325: !     Set up interpolations

327:       sh = 0.5*s
328:       th = 0.5*t
329:       sp = 0.5 + sh
330:       tp = 0.5 + th
331:       sm = 0.5 - sh
332:       tm = 0.5 - th
333:       shp(3,1) =   sm*tm
334:       shp(3,2) =   sp*tm
335:       shp(3,3) =   sp*tp
336:       shp(3,4) =   sm*tp

338: !     Set up natural coordinate functions (times 4)

340:       xo =  xl(1,1)-xl(1,2)+xl(1,3)-xl(1,4)
341:       xs = -xl(1,1)+xl(1,2)+xl(1,3)-xl(1,4) + xo*t
342:       xt = -xl(1,1)-xl(1,2)+xl(1,3)+xl(1,4) + xo*s
343:       yo =  xl(2,1)-xl(2,2)+xl(2,3)-xl(2,4)
344:       ys = -xl(2,1)+xl(2,2)+xl(2,3)-xl(2,4) + yo*t
345:       yt = -xl(2,1)-xl(2,2)+xl(2,3)+xl(2,4) + yo*s

347: !     Compute jacobian (times 16)

349:       xsj1 = xs*yt - xt*ys

351: !     Divide jacobian by 16 (multiply by .0625)

353:       xsj = 0.0625*xsj1
354:       if (xsj1.eq.0.0) then
355:          xsj1 = 1.0
356:       else
357:          xsj1 = 1.0/xsj1
358:       endif

360: !     Divide functions by jacobian

362:       xs  = (xs+xs)*xsj1
363:       xt  = (xt+xt)*xsj1
364:       ys  = (ys+ys)*xsj1
365:       yt  = (yt+yt)*xsj1

367: !     Multiply by interpolations

369:       ytm =  yt*tm
370:       ysm =  ys*sm
371:       ytp =  yt*tp
372:       ysp =  ys*sp
373:       xtm =  xt*tm
374:       xsm =  xs*sm
375:       xtp =  xt*tp
376:       xsp =  xs*sp

378: !     Compute shape functions

380:       shp(1,1) = - ytm+ysm
381:       shp(1,2) =   ytm+ysp
382:       shp(1,3) =   ytp-ysp
383:       shp(1,4) = - ytp-ysm
384:       shp(2,1) =   xtm-xsm
385:       shp(2,2) = - xtm-xsp
386:       shp(2,3) = - xtp+xsp
387:       shp(2,4) =   xtp+xsm

389:       end

391: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
392: !     int2d
393: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
394:       subroutine int2d(l,sg)
395: !-----[--.----+----.----+----.-----------------------------------------]
396: !     Purpose: Form Gauss points and weights for two dimensions

398: !     Inputs:
399: !     l       - Number of points/direction

401: !     Outputs:
402: !     sg(3,*) - Array of points and weights
403: !-----[--.----+----.----+----.-----------------------------------------]
404:       implicit  none
405:       PetscInt   l,i,lr(9),lz(9)
406:       PetscReal    g,third,sg(3,*)
407:       data      lr/-1,1,1,-1,0,1,0,-1,0/,lz/-1,-1,1,1,-1,0,1,0,0/
408:       data      third / 0.3333333333333333 /

410: !     2x2 integration
411:       g = sqrt(third)
412:       do i = 1,4
413:          sg(1,i) = g*lr(i)
414:          sg(2,i) = g*lz(i)
415:          sg(3,i) = 1.0
416:       end do

418:       end

420: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
421: !     ex54_psi - anusotropic material direction
422: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
423:       PetscReal function ex54_psi(x,y)
424:       implicit  none
425:       PetscReal x,y,theta
426:       common /ex54_theta/ theta
427:       ex54_psi = theta
428:       if (theta < 0.) then     ! circular
429:          if (y==0) then
430:             ex54_psi = 2.0*atan(1.0)
431:          else
432:             ex54_psi = atan(-x/y)
433:          endif
434:       endif
435:       end

437: !
438: !/*TEST
439: !
440: !   build:
441: !
442: !   test:
443: !      nsize: 4
444: !      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 -pc_gamg_esteig_ksp_max_it 5
445: !      requires: !single
446: !
447: !TEST*/