femwind/fortran/module_ndt_mult.f90

   1 module module_ndt_mult
   2
   3 use module_utils
   4
   5 contains
   6
   7 subroutine ndt_mult(                              &
   8     ifds, ifde, kfds,kfde, jfds, jfde,                       & ! fire grid dimensions
   9     ifms, ifme, kfms,kfme, jfms, jfme,            &
  10     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  11     ifts, ifte, kfts, kfte, jfts,jfte,             &
  12     kmat, lambda, y, r, siz, relres)
  13
  14 implicit none
  15
  16 !*** purpose: compute r = y - K*lambda and 2-norm of r
  17
  18 !*** arguments
  19
  20 integer, intent(in)::                             &
  21     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  22     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  23     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  24     ifts, ifte, kfts, kfte, jfts,jfte                            ! fire tile bounds
  25
  26 integer, parameter:: msize = 14
  27 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme,msize):: kmat  ! global stiffness matrix
  28 real,intent(in),  dimension(ifms:ifme,kfms:kfme,jfms:jfme):: lambda, y          ! input vectors
  29 real,intent(out), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: r          ! output vector
  30 real, intent(out):: siz, relres
  31 !*** local
  32
  33 integer:: i,j,k,ie,je,ke
  34 real:: t, res
  35
  36 !** executable
  37 ie = snode(ifte,ifde,+1)
  38 je = snode(jfte,jfde,+1)
  39 ke = snode(kfte,kfde,+1)
  40
  41 print *,'ndt_mult: loop upper bounds',ie,ke,je
  42
  43 siz = 0.
  44 res = 0.
  45
  46 do j=jfts,je
  47   do k=kfts,ke
  48     do i=ifts,ie
  49       r(i,k,j) = y(i,k,j) - ( &
  50         kmat(i-1,k-1,j-1,14)*lambda(i-1,k-1,j-1) +  &
  51         kmat(i  ,k-1,j-1,13)*lambda(i  ,k-1,j-1) +  &
  52         kmat(i+1,k-1,j-1,12)*lambda(i+1,k-1,j-1) +  &
  53         kmat(i-1,k-1,j  ,11)*lambda(i-1,k-1,j  ) +  &
  54         kmat(i  ,k-1,j  ,10)*lambda(i  ,k-1,j  ) +  &
  55         kmat(i+1,k-1,j  , 9)*lambda(i+1,k-1,j  ) +  &
  56         kmat(i-1,k-1,j+1, 8)*lambda(i-1,k-1,j+1) +  &
  57         kmat(i  ,k-1,j+1, 7)*lambda(i  ,k-1,j+1) +  &
  58         kmat(i+1,k-1,j+1, 6)*lambda(i+1,k-1,j+1) +  &
  59         kmat(i-1,k  ,j-1, 5)*lambda(i-1,k  ,j-1) +  &
  60         kmat(i  ,k  ,j-1, 4)*lambda(i  ,k  ,j-1) +  &
  61         kmat(i+1,k  ,j-1, 3)*lambda(i+1,k  ,j-1) +  &
  62         kmat(i-1,k  ,j  , 2)*lambda(i-1,k  ,j  ) +  &  ! K(I,J)*x(J)   I=(i-1,j,k)  storing upper triangle of K only, K(I,J) stored in another row
  63         kmat(i  ,k  ,j  , 1)*lambda(i  ,k  ,j  ) +  &  ! K(I,I)*x(I)   I=(i,j,k)
  64         kmat(i  ,k  ,j  , 2)*lambda(i+1,k  ,j  ) +  &  ! K(I,J)*x(J)   J=(i+1,j,k)
  65         kmat(i  ,k  ,j  , 3)*lambda(i-1,k  ,j+1) +  &  ! etc rest of the row I in upper triangle
  66         kmat(i  ,k  ,j  , 4)*lambda(i  ,k  ,j+1) +  &
  67         kmat(i  ,k  ,j  , 5)*lambda(i+1,k  ,j+1) +  &
  68         kmat(i  ,k  ,j  , 6)*lambda(i-1,k+1,j-1) +  &
  69         kmat(i  ,k  ,j  , 7)*lambda(i  ,k+1,j-1) +  &
  70         kmat(i  ,k  ,j  , 8)*lambda(i+1,k+1,j-1) +  &
  71         kmat(i  ,k  ,j  , 9)*lambda(i-1,k+1,j  ) +  &
  72         kmat(i  ,k  ,j  ,10)*lambda(i  ,k+1,j  ) +  &
  73         kmat(i  ,k  ,j  ,11)*lambda(i+1,k+1,j  ) +  &
  74         kmat(i  ,k  ,j  ,12)*lambda(i-1,k+1,j+1) +  &
  75         kmat(i  ,k  ,j  ,13)*lambda(i  ,k+1,j+1) +  &
  76         kmat(i  ,k  ,j  ,14)*lambda(i+1,k+1,j+1) )
  77       siz = max(siz,abs(y(i,k,j)))
  78       res = max(res,abs(r(i,k,j)))
  79     enddo
  80   enddo
  81 enddo
  82
  83 relres = res/max(tiny(siz),siz)
  84
  85
  86 end subroutine ndt_mult
  87
  88 end module module_ndt_mult
  89
  90