Actual source code: ex54f.F

petsc-3.5.2 2014-09-08
Report Typos and Errors
  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 <finclude/petscsys.h>
 23: #include <finclude/petscvec.h>
 24: #include <finclude/petscmat.h>
 25: #include <finclude/petscksp.h>
 26: #include <finclude/petscpc.h>
 27: #include <finclude/petscviewer.h>
 28: #include <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_CHARACTER,'-ne',ne,flg,ierr)
 59:       h = 2.d0/ne
 60:       M = (ne+1)*(ne+1)
 61:       theta = 90.d0
 62: !     theta is input in degrees
 63:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-theta',theta,flg,
 64:      $     ierr)
 65:       theta = theta / 57.2957795
 66:       eps = 1.d0
 67:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-epsilon',eps,flg,
 68:      $     ierr)
 69:       ki = 2
 70:       call PetscOptionsGetRealArray(PETSC_NULL_CHARACTER,'-blob_center',
 71:      $     blb,ki,flg,ierr)
 72:       if ( .not. flg ) then
 73:          blb(1) = 0.d0
 74:          blb(2) = 0.d0
 75:       else if ( ki .ne. 2 ) then
 76:          print *, 'error: ', ki,
 77:      $        ' arguments read for -blob_center.  Needs to be two.'
 78:       endif
 79:       call PetscOptionsGetBool(PETSC_NULL_CHARACTER,'-out_matlab',
 80:      $     out_matlab,flg,ierr)
 81:       if (.not.flg) out_matlab = PETSC_FALSE;

 83:       ev(1) = 1.d0
 84:       ev(2) = eps*ev(1)
 85: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 86: !     Compute the matrix and right-hand-side vector that define
 87: !     the linear system, Ax = b.
 88: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 89: !  Create matrix.  When using MatCreate(), the matrix format can
 90: !  be specified at runtime.
 91:       call MatCreate(PETSC_COMM_WORLD,Amat,ierr)
 92:       call MatSetSizes( Amat,PETSC_DECIDE, PETSC_DECIDE, M, M, ierr )
 93:       if ( npe == 1 ) then
 94:          call MatSetType( Amat, MATAIJ, ierr )
 95:       else
 96:          call MatSetType( Amat, MATMPIAIJ, ierr )
 97:       endif
 98:       call MatMPIAIJSetPreallocation(Amat,9,PETSC_NULL_INTEGER,6,
 99:      $     PETSC_NULL_INTEGER, ierr)
100:       call MatSetFromOptions( Amat, ierr )
101:       call MatSetUp( Amat, ierr )
102:       call MatGetOwnershipRange( Amat, Istart, Iend, ierr )
103: !  Create vectors.  Note that we form 1 vector from scratch and
104: !  then duplicate as needed.
105:       call MatGetVecs( Amat, PETSC_NULL_OBJECT, xvec, ierr )
106:       call VecSetFromOptions( xvec, ierr )
107:       call VecDuplicate( xvec, bvec, ierr )
108:       call VecDuplicate( xvec, uvec, ierr )
109: !  Assemble matrix.
110: !   - Note that MatSetValues() uses 0-based row and column numbers
111: !     in Fortran as well as in C (as set here in the array "col").
112:       thk = 1.d0              ! thickness
113:       nel = 4                   ! nodes per element (quad)
114:       ndf = 1
115:       call int2d(2,lint,sg)
116:       ix = 0
117:       do geq=Istart,Iend-1,1
118:          qj = geq/(ne+1); qi = mod(geq,(ne+1))
119:          x = h*qi - 1.d0; y = h*qj - 1.d0 ! lower left corner (-1,-1)
120:          if ( qi < ne .and. qj < ne ) then
121:             coord(1,1) = x;   coord(2,1) = y
122:             coord(1,2) = x+h; coord(2,2) = y
123:             coord(1,3) = x+h; coord(2,3) = y+h
124:             coord(1,4) = x;   coord(2,4) = y+h
125: ! form stiff
126:             ss = 0.d0
127:             do ll = 1,lint
128:                call shp2dquad(sg(1,ll),sg(2,ll),coord,shp,xsj,2)
129:                xsj = xsj*sg(3,ll)*thk
130:                call thfx2d(ev,coord,shp,dd,2,2,4,ex54_psi)
131:                j1 = 1
132:                do kj = 1,nel
133:                   a1 = (dd(1,1)*shp(1,kj) + dd(1,2)*shp(2,kj))*xsj
134:                   a2 = (dd(2,1)*shp(1,kj) + dd(2,2)*shp(2,kj))*xsj
135: c     Compute residual
136: c                  p(j1) = p(j1) - a1*gradt(1) - a2*gradt(2)
137: c     Compute tangent
138:                   i1 = 1
139:                   do ki = 1,nel
140:                      ss(i1,j1) = ss(i1,j1) + a1*shp(1,ki) + a2*shp(2,ki)
141:                      i1 = i1 + ndf
142:                   end do
143:                   j1 = j1 + ndf
144:                end do
145:             enddo

147:             idx(1) = geq; idx(2) = geq+1; idx(3) = geq+(ne+1)+1
148:             idx(4) = geq+(ne+1)
149:             if ( qj > 0 ) then
150:                call MatSetValues(Amat,4,idx,4,idx,ss,ADD_VALUES,ierr)
151:             else                !     a BC
152:                do ki=1,4,1
153:                   do kj=1,4,1
154:                      if (ki<3 .or. kj<3 ) then
155:                         if ( ki==kj ) then
156:                            ss(ki,kj) = .1d0*ss(ki,kj)
157:                         else
158:                            ss(ki,kj) = 0.d0
159:                         endif
160:                      endif
161:                   enddo
162:                enddo
163:                call MatSetValues(Amat,4,idx,4,idx,ss,ADD_VALUES,ierr)
164:             endif               ! BC
165:          endif                  ! add element
166:          if ( qj > 0 ) then      ! set rhs

168:             val = h*h*exp(-1.d2*((x+h/2)-blb(1))**2)*
169:      $           exp(-1.d2*((y+h/2)-blb(2))**2)
170:             call VecSetValues(bvec,1,geq,val,INSERT_VALUES,ierr)
171:          endif
172:       enddo
173:       call MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY,ierr)
174:       call MatAssemblyEnd(Amat,MAT_FINAL_ASSEMBLY,ierr)
175:       call VecAssemblyBegin(bvec,ierr)
176:       call VecAssemblyEnd(bvec,ierr)

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

182: !  Create linear solver context

184:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)

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

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

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

197:       call KSPSetFromOptions(ksp,ierr)

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

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


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

215:          call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Bvec',
216:      $        FILE_MODE_WRITE,viewer,ierr)
217:          call VecView(bvec,viewer,ierr)
218:          call PetscViewerDestroy(viewer,ierr)

220:          call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Xvec',
221:      $        FILE_MODE_WRITE,viewer,ierr)
222:          call VecView(xvec,viewer,ierr)
223:          call PetscViewerDestroy(viewer,ierr)

225:          call MatMult(Amat,xvec,uvec,ierr)
226:          val = -1.d0
227:          call VecAXPY(uvec,val,bvec,ierr)
228:          call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'Rvec',
229:      $        FILE_MODE_WRITE,viewer,ierr)
230:          call VecView(uvec,viewer,ierr)
231:          call PetscViewerDestroy(viewer,ierr)

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

259:       call VecDestroy(xvec,ierr)
260:       call VecDestroy(bvec,ierr)
261:       call VecDestroy(uvec,ierr)
262:       call MatDestroy(Amat,ierr)
263:       call KSPDestroy(ksp,ierr)
264:       call PetscFinalize(ierr)

266:       end

268: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: !     thfx2d - compute material tensor
270: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

272:       subroutine thfx2d(ev,xl,shp,dd,ndm,ndf,nel,dir)

274: c     Compute thermal gradient and flux

276:       implicit  none

278:       integer   ndm,ndf,nel,i
279:       PetscReal ev(2),xl(ndm,nel),shp(3,*),dir
280:       PetscReal xx,yy,psi,cs,sn,c2,s2,dd(2,2)

282: c      temp     = 0.0d0
283: c      gradt(1) = 0.0d0
284: c      gradt(2) = 0.0d0
285:       xx       = 0.0d0
286:       yy       = 0.0d0
287:       do i = 1,nel
288: c        gradt(1) = gradt(1) + shp(1,i)*ul(1,i)
289: c        gradt(2) = gradt(2) + shp(2,i)*ul(1,i)
290: c        temp     = temp     + shp(3,i)*ul(1,i)
291:         xx       = xx       + shp(3,i)*xl(1,i)
292:         yy       = yy       + shp(3,i)*xl(2,i)
293:       end do
294:       psi = dir(xx,yy)
295: c     Compute thermal flux
296:       cs  = cos(psi)
297:       sn  = sin(psi)
298:       c2  = cs*cs
299:       s2  = sn*sn
300:       cs  = cs*sn

302:       dd(1,1) = c2*ev(1) + s2*ev(2)
303:       dd(2,2) = s2*ev(1) + c2*ev(2)
304:       dd(1,2) = cs*(ev(1) - ev(2))
305:       dd(2,1) = dd(1,2)

307: c      flux(1) = -dd(1,1)*gradt(1) - dd(1,2)*gradt(2)
308: c      flux(2) = -dd(2,1)*gradt(1) - dd(2,2)*gradt(2)

310:       end

312: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313: !     shp2dquad - shape functions - compute derivatives w/r natural coords.
314: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
315:        subroutine shp2dquad(s,t,xl,shp,xsj,ndm)
316: c-----[--.----+----.----+----.-----------------------------------------]
317: c      Purpose: Shape function routine for 4-node isoparametric quads
318: c
319: c      Inputs:
320: c         s,t       - Natural coordinates of point
321: c         xl(ndm,*) - Nodal coordinates for element
322: c         ndm       - Spatial dimension of mesh

324: c      Outputs:
325: c         shp(3,*)  - Shape functions and derivatives at point
326: c                     shp(1,i) = dN_i/dx  or dN_i/dxi_1
327: c                     shp(2,i) = dN_i/dy  or dN_i/dxi_2
328: c                     shp(3,i) = N_i
329: c         xsj       - Jacobian determinant at point
330: c-----[--.----+----.----+----.-----------------------------------------]
331:       implicit  none
332:       integer   ndm
333:       real*8    xo,xs,xt, yo,ys,yt, xsm,xsp,xtm,xtp, ysm,ysp,ytm,ytp
334:       real*8    s,t, xsj,xsj1, sh,th,sp,tp,sm,tm, xl(ndm,4),shp(3,4)

336: c     Set up interpolations

338:       sh = 0.5d0*s
339:       th = 0.5d0*t
340:       sp = 0.5d0 + sh
341:       tp = 0.5d0 + th
342:       sm = 0.5d0 - sh
343:       tm = 0.5d0 - th
344:       shp(3,1) =   sm*tm
345:       shp(3,2) =   sp*tm
346:       shp(3,3) =   sp*tp
347:       shp(3,4) =   sm*tp

349: c     Set up natural coordinate functions (times 4)

351:       xo =  xl(1,1)-xl(1,2)+xl(1,3)-xl(1,4)
352:       xs = -xl(1,1)+xl(1,2)+xl(1,3)-xl(1,4) + xo*t
353:       xt = -xl(1,1)-xl(1,2)+xl(1,3)+xl(1,4) + xo*s
354:       yo =  xl(2,1)-xl(2,2)+xl(2,3)-xl(2,4)
355:       ys = -xl(2,1)+xl(2,2)+xl(2,3)-xl(2,4) + yo*t
356:       yt = -xl(2,1)-xl(2,2)+xl(2,3)+xl(2,4) + yo*s

358: c     Compute jacobian (times 16)

360:       xsj1 = xs*yt - xt*ys

362: c     Divide jacobian by 16 (multiply by .0625)

364:       xsj = 0.0625d0*xsj1
365:       if (xsj1.eq.0.0d0) then
366:          xsj1 = 1.0d0
367:       else
368:          xsj1 = 1.0d0/xsj1
369:       endif

371: c     Divide functions by jacobian

373:       xs  = (xs+xs)*xsj1
374:       xt  = (xt+xt)*xsj1
375:       ys  = (ys+ys)*xsj1
376:       yt  = (yt+yt)*xsj1

378: c     Multiply by interpolations

380:       ytm =  yt*tm
381:       ysm =  ys*sm
382:       ytp =  yt*tp
383:       ysp =  ys*sp
384:       xtm =  xt*tm
385:       xsm =  xs*sm
386:       xtp =  xt*tp
387:       xsp =  xs*sp

389: c     Compute shape functions

391:       shp(1,1) = - ytm+ysm
392:       shp(1,2) =   ytm+ysp
393:       shp(1,3) =   ytp-ysp
394:       shp(1,4) = - ytp-ysm
395:       shp(2,1) =   xtm-xsm
396:       shp(2,2) = - xtm-xsp
397:       shp(2,3) = - xtp+xsp
398:       shp(2,4) =   xtp+xsm

400:       end

402: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
403: !     int2d
404: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
405:       subroutine int2d(l,lint,sg)
406: c-----[--.----+----.----+----.-----------------------------------------]
407: c     Purpose: Form Gauss points and weights for two dimensions

409: c     Inputs:
410: c     l       - Number of points/direction

412: c     Outputs:
413: c     lint    - Total number of points
414: c     sg(3,*) - Array of points and weights
415: c-----[--.----+----.----+----.-----------------------------------------]
416:       implicit  none
417:       integer   l,i,lint,lr(9),lz(9)
418:       real*8    g,third,sg(3,*)
419:       data      lr/-1,1,1,-1,0,1,0,-1,0/,lz/-1,-1,1,1,-1,0,1,0,0/
420:       data      third / 0.3333333333333333d0 /

422: c     Set number of total points

424:       lint = l*l

426: c     2x2 integration
427:       g = sqrt(third)
428:       do i = 1,4
429:          sg(1,i) = g*lr(i)
430:          sg(2,i) = g*lz(i)
431:          sg(3,i) = 1.d0
432:       end do

434:       end

436: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
437: !     ex54_psi - anusotropic material direction
438: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
439:       PetscReal function ex54_psi(x,y)
440:       implicit  none
441:       PetscReal x,y,theta
442:       common /ex54_theta/ theta
443:       ex54_psi = theta
444:       if ( theta < 0. ) then     ! circular
445:          if (y==0) then
446:             ex54_psi = 2.d0*atan(1.d0)
447:          else
448:             ex54_psi = atan(-x/y)
449:          endif
450:       endif
451:       end