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: misk
439: !      args: -mat_coarsen_type misk -pc_gamg_aggressive_coarsening 0
440: !      output_file: output/ex54f_mis.out
441: !   test:
442: !      suffix: mis
443: !      args: -mat_coarsen_type mis
444: !      output_file: output/ex54f_mis.out
445: !
446: !TEST*/