Actual source code: ex54f.F90

petsc-3.11.2 2019-05-18
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*/


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

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

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

 85:       ev(1) = 1.0
 86:       ev(2) = eps*ev(1)
 87: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88: !     Compute the matrix and right-hand-side vector that define
 89: !     the linear system, Ax = b.
 90: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 91: !  Create matrix.  When using MatCreate(), the matrix format can
 92: !  be specified at runtime.
 93:       call MatCreate(PETSC_COMM_WORLD,Amat,ierr)
 94:       call MatSetSizes( Amat,PETSC_DECIDE, PETSC_DECIDE, M, M, ierr )
 95:       if ( size == 1 ) then
 96:          call MatSetType( Amat, MATAIJ, ierr )
 97:       else
 98:          call MatSetType( Amat, MATMPIAIJ, ierr )
 99:       endif
100:       call MatMPIAIJSetPreallocation(Amat,f9,PETSC_NULL_INTEGER,f6,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_VEC, 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.0              ! thickness
114:       nel = 4                   ! nodes per element (quad)
115:       ndf = 1
116:       call int2d(f2,sg)
117:       lint = 4
118:       ix = 0
119:       do geq=Istart,Iend-1,1
120:          qj = geq/(ne+1); qi = mod(geq,(ne+1))
121:          x = h*qi - 1.0; y = h*qj - 1.0 ! lower left corner (-1,-1)
122:          if ( qi < ne .and. qj < ne ) then
123:             coord(1,1) = x;   coord(2,1) = y
124:             coord(1,2) = x+h; coord(2,2) = y
125:             coord(1,3) = x+h; coord(2,3) = y+h
126:             coord(1,4) = x;   coord(2,4) = y+h
127: ! form stiff
128:             ss = 0.0
129:             do ll = 1,lint
130:                call shp2dquad(sg(1,ll),sg(2,ll),coord,shp,xsj,f2)
131:                xsj = xsj*sg(3,ll)*thk
132:                call thfx2d(ev,coord,shp,dd,f2,f2,f4,ex54_psi)
133:                j1 = 1
134:                do kj = 1,nel
135:                   a1 = (dd(1,1)*shp(1,kj) + dd(1,2)*shp(2,kj))*xsj
136:                   a2 = (dd(2,1)*shp(1,kj) + dd(2,2)*shp(2,kj))*xsj
137: !     Compute residual
138: !                  p(j1) = p(j1) - a1*gradt(1) - a2*gradt(2)
139: !     Compute tangent
140:                   i1 = 1
141:                   do ki = 1,nel
142:                      ss(i1,j1) = ss(i1,j1) + a1*shp(1,ki) + a2*shp(2,ki)
143:                      i1 = i1 + ndf
144:                   end do
145:                   j1 = j1 + ndf
146:                end do
147:             enddo

149:             idx(1) = geq; idx(2) = geq+1; idx(3) = geq+(ne+1)+1
150:             idx(4) = geq+(ne+1)
151:             if ( qj > 0 ) then
152:                call MatSetValues(Amat,f4,idx,f4,idx,ss,ADD_VALUES,ierr)
153:             else                !     a BC
154:                do ki=1,4,1
155:                   do kj=1,4,1
156:                      if (ki<3 .or. kj<3 ) then
157:                         if ( ki==kj ) then
158:                            ss(ki,kj) = .1*ss(ki,kj)
159:                         else
160:                            ss(ki,kj) = 0.0
161:                         endif
162:                      endif
163:                   enddo
164:                enddo
165:                call MatSetValues(Amat,f4,idx,f4,idx,ss,ADD_VALUES,ierr)
166:             endif               ! BC
167:          endif                  ! add element
168:          if ( qj > 0 ) then      ! set rhs
169:             val = h*h*exp(-100*((x+h/2)-blb(1))**2)*exp(-100*((y+h/2)-blb(2))**2)
170:             call VecSetValues(bvec,one,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)
204:       CHKERRA(ierr)


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

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

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

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

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

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

263:       end

265: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
266: !     thfx2d - compute material tensor
267: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268: !     Compute thermal gradient and flux

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

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

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

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

296: !      flux(1) = -dd(1,1)*gradt(1) - dd(1,2)*gradt(2)
297: !      flux(2) = -dd(2,1)*gradt(1) - dd(2,2)*gradt(2)

299:       end

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

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

327: !     Set up interpolations

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

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

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

349: !     Compute jacobian (times 16)

351:       xsj1 = xs*yt - xt*ys

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

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

362: !     Divide functions by jacobian

364:       xs  = (xs+xs)*xsj1
365:       xt  = (xt+xt)*xsj1
366:       ys  = (ys+ys)*xsj1
367:       yt  = (yt+yt)*xsj1

369: !     Multiply by interpolations

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

380: !     Compute shape functions

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

391:       end

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

400: !     Inputs:
401: !     l       - Number of points/direction

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

413: !     Set number of total points

415:       lint = l*l

417: !     2x2 integration
418:       g = sqrt(third)
419:       do i = 1,4
420:          sg(1,i) = g*lr(i)
421:          sg(2,i) = g*lz(i)
422:          sg(3,i) = 1.0
423:       end do

425:       end

427: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
428: !     ex54_psi - anusotropic material direction
429: !     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
430:       PetscReal function ex54_psi(x,y)
431:       implicit  none
432:       PetscReal x,y,theta
433:       common /ex54_theta/ theta
434:       ex54_psi = theta
435:       if ( theta < 0. ) then     ! circular
436:          if (y==0) then
437:             ex54_psi = 2.0*atan(1.0)
438:          else
439:             ex54_psi = atan(-x/y)
440:          endif
441:       endif
442:       end

444: !
445: !/*TEST
446: !
447: !   build:
448: !      requires: !libpgf90
449: !
450: !   test:
451: !      nsize: 4
452: !      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 -mg_levels_esteig_ksp_type cg
453: !      requires: !single
454: !
455: !TEST*/